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ABSTRACT 


The  study  of  the  formation  and  growth  of  thermal  plumes  is  motivated  by  the 
proposed  existence  of  such  plumes  in  the  Earth’s  mantle.  During  the  initial  stages  of 
plume  development,  a  plume  consists  of  a  large  buoyant  ball  trailing  a  narrow  feeder 
conduit.  This  study  presents  laboratory,  analytical,  and  numerical  models  of  this 
flow.  The  experimental  model  generates  the  plumes  using  a  heater  in  a  syrup  whose 
viscosity  is  highly  temperature-dependent.  The  resulting  data  provides  a  measure  of 
the  effectiveness  of  the  analytical  and  numerical  models.  The  analytical  model, 
based  on  mass  and  energy  conservation,  shows  a  significant  improvement  in  the  flow 
prediction  compared  to  previous  models.  The  numerical  model  uses  the  finite- 
element  method  to  produce  a  flow  solution  that  successfully  predicts  the  flow  to 
within  the  experimental  error. 


CHAPTER  1 


INTRODUCTION  AND  LITERATURE  REVIEW 

l.I  Introduction 

The  initiation  and  development  of  creeping  thermal  plumes  in  a  medium 
whose  viscosity  is  highly  temperature-dependent  presents  an  interesting  and 
challenging  problem  in  low  Reynolds  number  flow.  The  problem  gains  increased 
significance,  however,  when  considering  the  application  of  such  flow  to  geophysical 
fluid  dynamics  and  mantle  convection. 

Three  methods  of  investigation  were  used  to  researcn  this  problem: 
experimental,  numerical,  and  analytical.  A  brief  description  of  the  experimental 
apparatus  will  serve  to  provide  a  concrete  example  of  the  problem.  A  large  container 
was  filled  with  a  thick  corn  syrup  whose  viscosity  varies  greatly  with  temperature. 

A  small  heater  placed  in  the  bottom  of  the  container  provided  the  heat  input  to 
generate  a  thermally  buoyant  plume  which  formed  with  a  balloon-on-a-string  shape 
(figures  4c,  5c).  We  call  this  plume  structure  a  starting  plume.  The  goal  of  this 
project  was  to  predict  analytically  and  numerically  the  motion  of  the  fluid  in  the 
tank.  Of  particular  interest  was  the  rise  time  of  the  plume. 

This  problem  provides  a  first-approximation  model  of  deep  Earth-mantle 
plumes.  Morgan  (1971)  first  proposed  the  existence  of  mantle  plumes  extending  from 
the  deep  mantle  to  the  asthcnosphere.  While  several  regions  have  been  suggested  as 
the  source  of  mantle  plumes,  the  most  plausible  appears  to  be  the  D"  layer  (Stacey 
and  Loper,  1983).  Seismic  studies  have  identified  the  D"  layer  as  a  distinct  layer  at 
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the  base  of  the  mantle.  Yuen  and  Peltier  (1980)  used  linear  stability  analysis  to  show 
that  for  a  wide  range  of  conditions,  the  D"  layer  is  dynamically  unstable.  They 
postulated  that  this  instability  could  evolve  into  deep  mantle  plumes.  Stacey  and 
Loper  (1983)  modelled  the  D"  layer  as  a  dynamically  active  layer  with  a  flow  pattern 
that  is  distinct  from  that  in  the  lower  mantle.  This  layer,  maintained  by  the  heat 
flux  from  the  core,  provides  a  source  of  buoyant,  low  viscosity  fluid  for  plume 
activity.  At  the  high  pressures  and  temperatures  in  the  Earth’s  interior,  the  rock-like 
mantle  exhibits  fluid  properties.  A  key  characteristic  of  the  mantle  material  is  its 
highly  temperature-dependent  viscosity.  The  corn  syrup  in  the  experimental  tank, 
then,  represents  the  Earth’s  mantle.  The  heater  is  the  core  and  the  heated  fluid  layer 
just  above  the  heater  is  the  D"  layer.  While  this  crude  model  neglects  many 
complicating  factors  such  as  a  stably-stratified  mantle,  this  problem  serves  as  a  good 
initial  step  in  the  solution  of  the  larger,  more  complicated  problem. 

1.2  Literature  Review 

Most  of  the  previous  work  has  concentrated  on  specific  aspects  of  the  problem 
of  plume  initiation  and  growth;  few  have  attempted  to  integrate  these  aspects.  The 
efforts  may  be  grouped  into  two  categories:  experimental/analytical  investigations 
and  numerical  investigations. 

1.2.1  Experimental /Analytical  Investigations 

The  earliest  experimental/analytical  effort  appears  to  be  that  of  Whitehead 
and  Luther  (1975).  Many  of  the  results  were  summarized  in  Whitehead  (1986).  Their 
experiment  involved  injecting  a  less  dense,  less  viscous  fluid  into  a  denser,  more 
viscous  fluid.  The  viscosity  contrast  (Pi,r*e/l*fni»ii)  *^^0  fluids  was  6  x  10®. 

They  found  that  a  ball  formed  on  the  spout  of  the  injector,  and  then  lifted  off  from 
the  spout,  trailing  a  thin  feeder  conduit.  This  structure  is  very  similar  to  that 
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observed  in  the  research  experiment.  One  key  difference  in  the  experiments  is  that 
Whitehead  and  Luther  (1975)  used  two  immiscible  fluids,  so  no  diffusion  occurred. 
They  developed  an  analytical  model  for  speed,  height  and  radius  of  the  ball  of  the 
plume  as  a  function  of  time,  which  predicts  that  height  varies  as  t®^*.  Their 
experimental  results  showed  good  quantitative  agreement  with  the  model. 

A  similar  study  is  that  of  Olson  and  Singer  (1985).  They  also  continuously 
injected  a  less  dense,  less  viscous  fluid  into  a  more  viscous,  denser  fluid,  but  they 
formed  the  less  dense  fluid  by  diluting  the  denser  fluid  (glucose  syrup)  with  water 
forming  a  plume  of  miscible  fluid.  They  found  experimentally  that  height  varied  as 
t^^®  and  that  the  plume  rose  slower  than  for  the  Whitehead  and  Luther  (1975)  case. 
This  slower  rate  may  be  due  to  diffusion.  Note  that  the  compositional  diffusion  of 
the  Olson  and  Singer  (1985)  experiment  is  very  slow  relative  to  the  thermal  diffusion 
of  the  experiment  in  the  current  study. 

Diffusion  of  buoyancy  is  a  key  factor  in  the  thermal  plume  problem,  as 
indicated  in  the  research  of  Morris  (1982)  and  Ansari  and  Morris  (1985).  These 
studies  isolated  the  effects  of  strongly  temperature-dependent  viscosity  on  slow  flow 
past  a  rigid  sphere.  These  investigations  were  motivated  by  the  observation  that  a 
sphere  of  hot  magma  would  require  a  time  equivalent  to  the  Earth’s  age  to  move 
through  the  lithosphere  in  a  Stokes-flow  regime.  Morris  showed  that  if  a  rigid 
sphere  is  maintained  at  a  high  temperature  in  a  medium  with  highly  temperature- 
dependent  viscosity,  then  it  will  transfer  heat  to  the  surrounding  medium  and  create 
a  lower  viscosity  path  for  the  sphere  to  follow.  In  this  ’lubrication  limit’,  the  sphere 
moves  faster  than  in  Stokes  flow.  In  fact,  Morris  found  for  the  lubrication  limit  that 
a  smaller  sphere  rises  faster  than  a  larger  sphere,  an  effect  opposite  to  Stokes  flow. 
This  result  was  confirmed  experimentally  by  Ribe  (1983).  Ribe  also  determined  that 
the  ascent  speed  will  increase  by  25%  (for  Stokes  flow)  to  100%  (lubrication  limit) 
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for  each  order-of -magnitude  change  in  viscosity  contrast.  These  results  show  that  a 
sphere  of  magma  would  rise  faster  than  predicted  by  Stokes  flow,  if  the  temperature 
of  the  sphere  were  constant.  In  reality,  such  a  sphere  would  cool  quickly  and  stop  its 
upward  motion. 

Griffiths  (1986a-d)  made  an  experimental  investigation  that  also  concentrated 
on  the  effects  of  diffusion.  In  contrast  to  the  continuously  injected  plumes  of 
Whitehead  and  Luther  (1975)  and  Olson  and  Singer  (1985),  he  injected  a  fixed 
volume  of  warmer  oil  into  a  larger  container  of  the  same  oil  at  a  lower  temperature. 
These  ’thermals’  grew  in  size  due  to  entrainment  of  the  ambient  fluid,  while  the 
speed  decreased  with  time.  By  assuming  that  all  of  the  heated  ambient  fluid  is 
entrained  (i.e.,  buoyancy  is  conserved),  Griffiths  developed  analytical  expressions 
indicating  that  height  varied  as  t^^^  and  speed  varied  as  This  rise  law  agreed 
qualitatively  and  quantitatively  with  his  experimental  results. 

Recent  work  by  Griffiths  and  Campbell  (1990)  more  closely  models  the 
proposed  experimental  study.  This  work  also  represents  the  first  attempt  to  integrate 
ail  aspects  of  the  thermal-plume  problem.  They  injected  a  continuous  stream  of 
warmer,  less  viscous  fluid  into  an  overlying  layer  of  more  dense,  more  viscous  fluid. 
By  assuming  that  all  of  the  heated  ambient  fluid  is  entrained,  they  developed  a  rise 
law  predicting  that  the  height  of  the  plume  would  vary  as  t®/^  Although  this 
assumption  is  not  physically  realistic,  they  stated  that  this  rise  law  adequately 
predicted  their  experimental  results,  but  presented  no  data  to  support  this  statement. 
1.2.2  Numerical  Investigations 

Just  as  the  previous  experimental/analytical  work  isolated  various  factors  of 
the  plume  flow,  much  of  the  numerical  work  emphasized  selected  physical 
phenomena.  Some  of  the  numerical  work,  however,  attempts  to  attack  the  full 
problem  in  methods  similar  to  the  proposed  research.  Table  1  summarizes  the 
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Table  1.  Comparison  of  related  numerical  studies. 


Name 

Scheme 

Domain 

Ra 

Pc/Ph 

Note 

Parmentier  (75) 

FDM 

axisym 

10® 

10*- 

10® 

steady 

state 

Daly  and 

Raefsky  (85) 

FEM 

axisym 

0 

10® 

Schubert  and 
Anderson  (85) 

FEM 

planar 

107 

n/a 

constant 

viscosity 

Jarvis (84) 

FDM 

planar 

10® 

n/a 

constant 

viscosity 

Boss  and 

Sacks (85) 

FDM 

planar 

10^ 

20 

stream 

funct¬ 

ion 

Olson,  Schubert, 
Anderson  (87) 

FEM 

planar 

10^ 

1 

o  o 

Zhao  and 

Yuen  (87) 

FEM 

planar 

axisym 

10^ 

100 

viscous 

heating 

characteristics  of  previous  related  numerical  work.  Both  finite-difference  methods 
and  finite-element  studies  are  described  in  the  literature.  In  general,  the  more  recent 
works  prefer  the  finite-element  method  due  to  the  grid  flexibility  which  allows  high 
resolution  in  the  boundary  layers  that  commonly  occur  in  such  flows.  Sato  and 
Thompson  (1976)  discuss  some  of  the  specific  advantages  of  the  finite-clement 
method  for  this  application.  Parmentier,  et  al.  (1975)  made  the  first  significant 
numerical  study  on  the  structure  of  mantle  plumes.  They  used  a  finite-difference 
method,  which  incorporated  the  Boussinesq  approximation  and  an  infinite-Prandtl- 
number  (Pr)  assumption,  to  solve  the  steady-state  axisymmctric  problem.  They 
solved  the  energy  and  vorticity  transport  equations  containing  two  parameters, 
Rayleigh  number  (Ra)  and  viscosity  contrast  (p^)  where: 


where: 
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V^K  K 

a  is  the  coefficient  of  thermal  expansion 

g  is  the  acceleration  of  gravity 

AT  is  the  temperature  difference 

h  is  the  height 

K  is  the  thermal  diffusivity 

Pf  is  the  reference  dynamic  viscosity  of  the  cold 

fluid 

p  is  the  dynamic  viscosity  of  the  hot  fluid 
is  the  reference  kinematic  viscosity  of  cold 
fluid. 

Although  they  assumed  flow  was  restricted  to  the  upper  mantle  (700  km),  some  of 
their  conclusions  are  equally  valid  for  deep  mantle  plumes.  They  found,  for 
instance,  that  the  medium  must  be  heated  at  the  base  for  plumes  to  occur;  internal 
heating  would  not  produce  the  narrow  structure.  They  also  showed  that  variable 
viscosity  was  necessary  for  narrow  plumes.  Finally,  they  investigated  the  effect  of 
varying  the  Rayleigh  number  from  10^  to  10®.  They  showed  that  as  the  Rayleigh 
number  increased,  the  plume  neck  width  decreased. 

Daly  and  Raefsky  (1985)  performed  a  numerical  study  of  the  situation 
described  asymptotically  by  Morris  (1982)  and  experimentally  by  Ribe  (1983),  i.e., 
the  motion  of  a  rigid  sphere  in  a  highly  temperature-dependent  medium.  They  used 
a  finite-element  method  on  the  momentum  and  energy  equations,  incorporating  some 
of  the  same  assumptions  as  Parmentier,  et  al.  (1975),  specifically  the  Boussinesq 
approximation  and  infinite-Prandtl-number  assumption,  but  they  included  time 
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dependency  in  the  energy  equation.  They  also,  incidentally,  used  a  primitive- 
variable  approach  with  a  penalty-function  formulation  in  place  of  the  vorticity 
approach.  Their  results  showed  excellent  agreement  with  the  results  of  Morris 
(1982),  Ribe  (1983),  and  Ansari  and  Morris  (1983).  Additionally,  they  were  able  to 
investigate  a  wider  parameter  range:  viscosity  variation  up  to  10®  and  Peclet  number 
(Pe)  from  10"^  to  10®,  where 


K 


where: 

U  is  the  speed  of  the  sphere 
a  is  the  sphere  radius. 

Daly  and  Raefsky  (1985)  assumed  the  fluid  was  mechanically  driven  for  their  work. 
They  concluded  that  although  the  drag  reduces  with  viscosity  variation  (as  compared 
to  Stokes  flow  with  constant  viscosity),  the  heat  transport  is  also  more  efficient,  so  a 
blob  of  magma  will  not  rise  farther  than  a  diameter  or  two  through  the  mantle  or 
lithosphere  without  cooling  to  ambient  levels. 

Many  studies  concentrate  on  overall  mantle  convection  but  don’t  investigate 
plume  structure  in  detail.  Several  studies,  while  concentrating  on  overall  mantle 
convection,  mention  the  role  of  plumes  in  their  results.  Schubert  and  Anderson 
(1985)  investigated  two-dimensional  convection  with  heating  from  below  as  well  as 
internal  heating  at  high  Rayleigh  numbers  using  the  finite-element  method.  They 
compare  their  results  favorably  with  earlier  finite-difference  studies  of  Jarvis  and 
Peltier  (1982).  While  they  made  progress  in  using  a  high  Rayleigh  number  that 
would  be  typical  of  whole-mantle  convection,  their  model  assumed  constant  viscosity 
in  the  mantle.  Jarvis  (1984)  was  able  to  study  flow  at  even  higher  Rayleigh  numbers 
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using  a  finite-difference  method  and  a  stream-function  approach  to  the  two- 
dimensional  problem,  but  he  also  assumed  constant  viscosity.  He  also  concentrated 
on  behavior  of  developed  plumes,  not  plume  initiation. 

Four  papers  are  most  closely  related  to  the  numerical  approach  used  in  this 
research.  They  each  make  the  Boussinesq  approximation  and  infinite  Prandtl 
number  assumption  in  a  two-dimensional  planar  domain.  As  demonstrated  by  the 
common  occurrence  of  these  assumptions  in  the  papers  described  earlier,  these 
assumptions  are  well  accepted  for  models  of  mantle  convection.  Boss  and  Sacks 
(1985)  solve  a  stream-function  formulation  of  the  problem  with  a  finite-difference 
method.  They  initiated  plumes  by  introducing  a  perturbation  into  broad-scale 
mantle  convection  as  well  as  a  static  mantle,  in  an  attempt  to  determine  which  core¬ 
mantle  boundary  perturbations  produced  plumes.  They  concluded  that  deep-mantle 
plumes  require  50-100  my  to  penetrate  the  mantle.  This  time  scale  was  confirmed  by 
Christensen’s  (1984)  study.  Olson,  Schubert  and  Anderson  (1987)  use  a  finite-clement 
code  to  model  plume  initiation  in  the  lower  mantle,  with  viscosity  contrasts  of  lO’ 
and  10^  Zhao  and  Yuen  (1987)  introduced  another  phenomenon  by  studying  the 
effects  of  adiabatic  and  viscous  heating  on  plumes.  They  also  are  the  only  authors  to 
use  an  axisymmetric  finite-element  formulation  in  addition  to  a  two-dimensional 
Cartesian  formulation  for  the  time-dependent  problem.  They  found  that  for  a 
uniform  thermal  expansivity,  a,  an  increase  in  the  adiabatic  and  viscous  heating 
contributions  impedes  plume  growth.  Only  with  a  decreasing  with  depth  docs  plume 
structure  return  in  this  model.  They  also  observed  that  the  same  plume  structure 
appears  in  the  Cartesian  two-dimensional  and  axisymmetric  cases,  the  only 
quantitative  difference  being  that  the  temperature  in  the  center  of  the  plume  is 
slightly  higher  in  the  two-dimensional  case  than  in  the  axisymmetric  case. 


CHAPTER  2 


EXPERIMENTAL  INVESTIGATION 

2.1  Purpose 

The  goal  of  the  experimental  study  was  to  measure  the  rise  speed  and  to 
observe  the  development  of  starting  thermal  plumes.  As  indicated  in  the  literature 
review,  a  key  aspect  of  this  study  was  the  combination  of  the  buoyancy  effects  seen 
in  the  continuously-fed  compositional  plume  studies  of  Whitehead  and  Luther  (1975) 
and  Olson  and  Singer  (1985)  with  the  entrainment  effects  demonstrated  by  the  pulse- 
fed  ’thermals’  of  Griffiths  (1986a-d).  Griffiths  and  Campbell  (1990)  first  studied 
this  combination  of  buoyancy  and  diffusive  effects  experimentally.  The  current 
experimental  work  differs  from  theirs  in  two  ways:  the  plumes  are  generated  via 
thermal  input  rather  than  mass  input  and  the  viscosity  contrast  between  the  fluids  is 
larger  than  their  case.  This  work  at  higher  viscosity  contrasts  may  be  helpful  in 
drawing  conclusions  about  plume  behavior  at  the  very  large  viscosity  contrasts 
which  may  occur  in  the  mantle. 


2.2  Apparatus 

The  spherical  caps  of  thermally  generated  starting  plumes  have  larger 
diameter  than  those  of  compositionally  generated  plumes  (Olson  and  Singer,  1985). 
This  difference  is  due  to  the  larger  value  of  thermal  diffusivity  as  compared  to 
compositional  diffusivity  and  the  resulting  entrainment  described  by  Griffiths 
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(1986a-d).  Thus,  compositional  plumes  may  be  studied  in  a  smaller  container  with 
minimal  wall  effects.  For  the  present  study,  wall  effects  are  minimized  for  thermal 
plumes  by  using  a  large  tank  (figure  1). 

The  tank  (designed  to  be  the  maximum  size  that  will  move  through  a  standard 
doorway)  has  dimensions  78.8  X  78.8  X  69.4  cm  high.  The  plume  was  generated  by  a 
62.2  mm  diameter  thermo-foil  heater  which  has  265  ohm  resistance.  The  heater  was 
mounted  on  a  metal  block;  a  7  mm  thick  layer  of  plexiglass  separated  the  heater 
from  the  block.  Vertical  length  scales  on  the  front  and  rear  of  the  tank  provided  a 
means  of  measuring  the  plume  height  from  the  photographic  record  of  the 
experiment. 

Table  2.  Physical  properties  of  syrup. 


thermal  diffusivity 

2.2  x  10‘*  cm*/scc 

specific  heat 

.7  cal/gm/deg 

thermal  expansion  (S)22C 

3.5  X  10-“*  deg-^ 

thermal  expansion  (®40C 

4.5  X  10'^  deg*^ 

density  (5)-26.1C 

1.45  g/cm’ 

density  @.1C 

1.43  g/cm® 

density  @25C 

1.42  g/cm’ 

density  @100C 

1.39  g/cm® 

dynamic  viscosity  (a>-26.1C 

2.59  X  10'°  mPa-sec 

dynamic  viscosity  (S».1C 

7.4  X  10®  mPa-sec 

dynamic  viscosity  (®25C 

91 104  mPa-sec 

dynamic  viscosity  @100C 

223.9  mPa-sec 

The  fluid  used  in  the  tank  was  ADM  36/43  industrial  corn  sweetener.  A 
summary  of  some  of  its  important  physical  properties  is  given  in  table  2.  This  syrup 
has  a  highly  temperature-dependent  viscosity  (figure  2),  a  key  characteristic  of 


Figure  1.  Experimental  apparatus. 
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mantle  material.  The  table  also  indicates  a  significant  change  in  the  thermal 
expansion  coefficient  with  temperature.  For  simplification  in  analytical  and 
numerical  modelling,  however,  we  assumed  a  constant  value  for  each  case.  Another 
useful  property  of  this  syrup  is  that  when  heated  in  the  experiment,  many  small 
vapor  bubbles  form,  apparently  the  residual  of  larger  bubbles  that  formed  on  the 
surface  of  the  heater  and  quickly  collapsed  as  they  migrated  away  from  the  heater. 
These  remnants  serve  as  excellent  tracers  for  flow  visualization,  yet  they  do  not 
appear  to  affect  the  flow  properties.  This  observation  has  also  been  noted  by  Ansari 
and  Morris  (1985).  If  these  bubbles  occupied  a  large  percentage  volume  of  the 
plume,  then  they  would  significantly  reduce  the  density  in  the  plume.  The 
experimental  photographs  indicate,  however,  that  these  bubbles  probably  do  not 
occupy  a  significant  volume  percentage.  This  assumption  is  based  on  the  observation 
that  one  can  see  through  the  syrup/bubble  mixture  fairly  easily.  A  large  percentage 
of  bubbles  would  obscure  the  view.  As  a  result,  we  do  not  account  for  bubble 
density  in  the  plume  density  calculations. 

2.3  Procedure 

The  continuous-feed  experiments  involved  setting  the  heater  to  full  power  in 
the  tank  of  quiescent,  isothermal  syrup.  The  record  of  the  ball  shape,  volume  and 
height  as  a  function  of  time  was  kept  photographically  (figures  4c,  5c,  and  6b).  The 
experiments  were  run  repeatedly  for  various  viscosity  contrasts,  which  were  induced 
by  changing  the  ambient  temperature  of  the  tank.  Three  convenient  environments 
existed  to  run  the  experiment  with  different  viscosity  contrasts.  One,  of  course,  was 
room  temperature.  This  environment  had  an  ambient  temperature  of  24.8C  and  this 
temperature  generated  a  viscosity  contrast  (Pcoid/Phot)  406.9.  The  other  two 
environments  come  as  a  result  of  access  to  the  National  Science  Foundation’s 


log(viscosity)  log(mPa-sec) 


Figure  2.  Viscosity  variation  with  temperature  for  syrup.  Least  squares  fit 
equation  is  logio(^)-20.6-14.99(1000/T)+3.0762(1000/T)* 
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Antarctic  Marine  Research  Facility  on  the  campus  of  Florida  State  University.  This 
facility  has  large  cold  storage  areas  that  will  accommodate  the  tank.  One  storage 
area  is  maintained  at  .1C,  generating  a  viscosity  contrast  of  3.3  X  10^.  A  smaller 
storage  area,  maintained  at  -26. 1C,  generates  a  viscosity  contrast  of  1.1  X  10®.  For 
comparison,  note  that  the  Griffiths  and  Campbell  (1990)  experiments  used  a 
maximum  viscosity  contrast  of  333. 

In  addition  to  the  continuous-feed  cases  that  were  run  at  different  viscosities, 
a  control  pulse-feed  case  was  run  at  ambient  temperature  of  .1C.  In  this  case,  the 
heater  was  turned  off  as  soon  as  the  sphere  formed  above  the  heater,  just  before  the 
time  when  the  ball  would  rise  away  from  the  heater,  trailing  a  narrow  neck  in  the 
continuous-feed  case.  (We  define  this  time  as  the  liftoff  time.)  This  case  was  very 
similar  to  the  experimental  study  of  ’thermals’  made  by  Griffiths  (I986a-d)  in  which 
he  injected  a  fixed  volume  of  warmer  fluid  into  the  ambient  fluid.  Comparison  of 
the  results  to  those  of  Griffiths  (1986a-d)  provided  a  good  test  of  the  experimental 
apparatus  and  procedure. 


2.4  Results 

The  results  of  the  control  pulse-feed  case  are  presented  in  figure  3a.  A 
photograph  of  the  typical  structure  of  the  thermal  is  shown  in  figure  3b.  The 
general  agreement  of  these  results  with  those  of  Griffiths  (1986a-d)  indicates  that 
the  experimental  apparatus  and  procedure  were  adequate.  Furthermore,  this 
agreement  suggests  that  plume  formation  via  thermal  input  (heater)  occurs  similarly 
to  plume  formation  via  injection  of  warmer  fluid  as  in  Griffiths  (1986a-d)  and 
Griffiths  and  Campbell  (1990). 

Typical  results  for  the  continuous-feed  cases  for  all  three  environments  are 
given  in  figures  4,  5,  and  6.  The  most  accurate  measurements  were  those  of  height 
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versus  time.  These  results  are  given  in  figures  4a,  5a,  and  6a.  The  measurements  of 
volume  versus  time  are  much  less  accurate;  these  results  are  given  in  figures  4b  and 
5b.  Corresponding  photographs  of  the  typical  structures  in  each  environment  are 
shown  in  figures  4c,  5c,  and  6b.  Along  with  the  experimental  results,  plots  of  the 
rise  laws  of  Whitehead  and  Luther  (1975),  Olson  and  Singer  (1985),  and  Griffiths  and 
Campbell  (1990)  are  included  in  figures  4a,  4b,  5a,  5b,  and  6a.  Note  the  error  bars  on 
the  graphs.  The  agreement  of  these  rise  laws  with  the  experimental  data  depends 
largely  upon  where  the  plot  begins.  In  order  to  have  a  fair  comparison,  the  plot  of 
each  rise  law  starts  at  the  liftoff  time  and  passes  through  the  first  two  observational 
points.  These  plots  provide  a  graphical  means  of  comparing  the  effectiveness  of  the 
rise  laws.  Additionally,  table  3  provides  a  quantitative  means  of  comparing  these 
models  (using  height  vs.  time  data).  By  computing  the  experimental  standard 
deviation  of  the  least-squares  fit  for  each  rise  law,  a  numerical  value  may  be 
assigned  which  indicates  how  well  each  rise  law  models  the  experimental  result. 

Note  that  these  results  represent  curves  that  are  different  from  those  in  figures  4a, 
5a,  and  6a.  The  continuous-feed  case  with  viscosity  contrast  of  1.1  X  10®  shown  in 
figure  6a  and  6b  displayed  a  different  morphology  than  previous  cases;  the  ball 
never  lifted  away  from  the  heater. 

The  graphical  comparison  of  these  plots  with  the  experimental  data  and  the 
results  presented  in  table  3  lead  to  a  few  observations.  First,  the  rise  law  of 
Whitehead  and  Luther  (1975)  is  the  poorest  model  of  the  experimental  data.  This 
poor  performance  is  expected,  however,  since  Whitehead  and  Luther  did  not 
incorporate  any  diffusive  effects  in  their  model,  nor  should  they  since  they  were 
dealing  with  immiscible  fluids.  Griffiths’  (1986a-d)  work  clearly  shows  that 
diffusion  and  subsequent  entrainment  of  ambient  fluid  accelerate  the  growth  of  the 
sphere  and  slow  the  sphere  rise  when  compared  to  the  no  diffusion  case,  so 
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Figure3a.  Results  of  pulse  feed  experiment. 
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Figure  3b.  Typical  structure  of  pulse  feed  diapir. 
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Figure  4a.  Height  vs.  time,  T,^jj«25C. 
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rc  4c.  Typical  structure  for  Tjj^j^=25C  ease. 
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Figure  5b.  Volume  vs.  time,  T 
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Figure  5c.  Typical  structure  of  T 
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igure  6a.  Height  vs.  time,  T^^«-26.1C. 
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r  the  thermal  plume  will  result  in  larger  rise  speeds.  The 
75)  law,  consequently,  shows  an  increasing  speed  since  the 
Dy  the  slowing  effects  of  diffusion.  The  volume  of  the  ball 
mailer  for  this  case  since  their  model  only  accounts  for  the 
nduit,  i.  e.,  no  ambient  material  is  entrained.  The  rise  law 
)  performs  better  than  that  of  Whitehead  and  Luther 
iequate.  Their  rise  law  is  based  on  experimental 
}nal  analysis,  and  does  not  consider  the  dynamics  of  the 
plumes  were  generated  using  a  compositional  difference  of 
did  not  include  compositional  diffusion  is  a  factor  in  their 
eir  experimental  results,  howeve*-,  shewed  some  slowing  and 
to  the  Whitehead  and  Luthei  (1975)  rise  law.  These 
small  amount  of  compositional  diffusion  between  the  two 
lal  diffusion  is  a  slow  process  relative  to  thermal  diffusion, 
ly,  the  best  model  to  date  is  that  of  Griffiths  and  Campbell 
5  the  effect  of  thermal  diffusion.  A  key  assumption  of 
hat  all  of  the  heat  diffusing  from  the  bail  is  re-entrained  in 
^es  the  ball.  Their  model  predicts  a  higher  speed  and  a 
time  than  the  experimental  results  indicate,  as  shown  in 

ble  3  lists  the  minimum  standard  deviation  for  an 
exponential  rise  law  is  one  where  height  -  time  Each  of 
ws  is  an  exponential  rise  law  »  5/3  for 
5  for  Olson/Singer,  and  {  =  5/4  for  Griffiths/Campbell), 
the  exponential  rise  law  that  predicts  the  experimental 
as  i  -  1.047.  For  the  T,„b  *  O.IC  case,  5  -  1.125,  and  for 
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Table  3.  Quantitative  model  comparison. 


MODEL 

STANDARD  DEVIATION 

T.„,b=25C 

T.„b=0.1C 

Whitehead/Luther 

1.44602 

1.02836 

3.47298 

Olson/Singer 

.94595 

.60531 

2.87026 

Griffiths/Campbell 

.69916 

.40745 

2.5039 

optimal  exponential 

.5304 

.3348 

.17538 

the  =  -26.1C  case,  C  =  .3743.  The  Olson  and  Singer  (1985)  experimental  rise  law 
gives  a  54%  improvement  (for  the  T^j,  =  25c  case)  and  61%  improvement  (for  the 
Tjmb  =  case)  over  the  Whitehead  and  Luther  (1975)  rise  law  when  compared  to 
this  optimal  exponential  case.  The  Griffiths  and  Campbell  (1990)  addition  of 
diffusion  to  the  model  shows  a  larger  improvement  in  the  model’s  performance  when 
compared  to  these  optimal  exponential  rise  laws.  Their  model  accounted  for  81.6% 
(in  the  =  25C  case)  and  89.5%  (in  the  =  O.IC  case)  of  the  difference 
between  the  Whitehead  and  Luther  (1975)  model  and  the  best  exponential  rise  law. 
Note  that  for  the  =  -26.1C  case,  the  Olson  and  Singer  (1985)  rise  law  shows  an 
18.2%  improvement  over  the  Whitehead  and  Luther  (1975)  prediction.  The  Griffiths 
and  Campbell  (1990)  model  gives  an  improvement  of  29%,  but  the  relatively  large 
standard  deviations  in  all  the  models  indicate  that  none  of  the  previously  developed 
rise  laws  is  adequate  for  this  case. 

These  results  show  the  importance  of  including  diffusive  effects  in  an 
analytical  model  of  plume  rise.  They  also  indicate  the  possibility  of  improving  the 
model  by  relaxing  some  of  the  assumptions  used  by  Griffiths  and  Campbell  (1990). 
Figure  6b  suggests  that  their  assumption  of  total  heat  entrainment  may  be  too  strong. 
In  this  T,^b  *  *26. 1C  case,  clearly  much  of  the  heat  input  to  the  plume  is  lost  from 
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the  ’ball’.  This  result  suggests  that,  to  a  lesser  degree,  some  heat  may  be  lost  from  the 
ball  in  the  other  two  cases.  A  logical  step  for  improving  their  model  is  to  relax  this 
assumption. 


CHAPTER  3 


ANALYTICAL  MODELLING 

3.1  Key  Physical  Interactions 

One  of  the  major  reasons  for  performing  an  experimental  inyestigation  of 
this  problem  was  to  deyclop  some  sense  of  the  important  physical  interactions  that 
occur  in  plume  initiation  and  deyelopment.  The  experimental  results  clearly 
indicate  that  one  important  interaction  is  thermal  diffusion  and  entrainment.  The 
hot  fluid  rising  from  the  heater  heats  the  ambient  fluid  and  entrains  some  of  it  into 
the  sphere.  Griffiths  (1986a-d)  showed  that  this  entrainment  alone  causes  the  sphere 
to  grow  in  size  and  slow  in  speed  as  time  progresses.  In  the  experimental  results  of 
chapter  2  (figures  4a  and  5a),  the  performance  of  the  model  that  included  diffusion 
(Griffiths  and  Campbell,  1990)  was  superior  to  the  models  that  neglected  diffusion 
(Whitehead  and  Luther,  1975;  Olson  and  Singer,  1985).  Another  important 
interaction  that  the  experiment  illustrated  was  heat  loss.  In  the  T^j,  =  -26.1C  case 
(figure  6b),  a  large  fraction  of  heat  was  lost  from  the  ’ball’.  This  factor  may  also 
play  a  significant  role  in  the  cases  with  smaller  viscosity  contrast. 

The  complications  introduced  by  these  interactions  make  analytical 
investigation  difficult.  For  example,  boundary  layer  analysis  is  complicated  by  the 
fact  that  the  thermal  layer  is  wound  into  the  interior  of  the  ball.  While  these 
interactions  make  many  types  of  analysis  difficult  or  impossible,  they  also  force  a 
very  fundamental  approach  to  the  study.  The  goal  of  this  analysis  is  to  produce  a 
first  order  or  ’engineering’  approximation  to  the  fluid  behavior.  This  approximation 
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will  be  tested  against  the  experimental  results  presented  in  chapter  two. 

3.2  System  Description 

Figure  7  presents  a  sketch  used  to  define  the  system  under  study.  The  system 
consists  of  that  portion  of  the  (radius  c)  ball  of  the  plume  marked  by  the  stippled 
region.  The  center  of  the  ball  serves  as  the  origin  of  the  frame  of  reference  for  the 
system,  with  the  far-ficld  velocity,  U,  approaching  the  bail  from  above.  A  dividing 
streamline,  marked  with  the  letter  D,  separates  the  ambient  fluid  that  is  entrained 
from  that  which  is  left  behind.  The  system  is  defined  so  that  all  entrained  fluid 
enters  the  system  through  the  annular  region  marked  by  B.  The  warm  fluid  from  the 
heater  enters  the  system  through  the  annular  region  marked  by  A.  The  symbol  C 
represents  the  entrained-layer  thickness,  4  represents  the  thermal-boundary-layer 
thickness,  F  is  the  (constant)  source  volume  flux  into  the  system,  r^  is  the  radius  of 
the  feeder  conduit,  and  T,  is  the  temperature  of  the  source  material.  T^i,,  and 
Pe  represent  the  (constant)  quantities  of  ambient  temperature,  density,  and  viscosity. 

3.3  Assumptions 

With  the  system  so  defined,  make  the  following  simplifying  assumptions  for 
the  analytical  model: 

(1)  Heat  loss  in  the  conduit  is  negligible  compared  to  heat  loss  in  the 
ball. 

(2)  The  temperature,  density,  and  viscosity  inside  the  ball  may  be 
approximated  by  average  values.  Let  T,  p,  and  p  represent  the  average 
temperature,  density,  and  viscosity,  respectively,  of  the  fluid  inside 
the  system. 

(3)  The  temperature  in  the  thermal  layer  of  thickness  4  varies  only  in 


Figure  7.  System  used  in  analytical  solution 


32 


the  radial  direction.  Furthermore,  the  temperature  profile  across  this 
layer  is  linear. 

(4)  Negligible  differences  exist  between  the  values  of  density  and 
specific  heat  in  the  heated  and  ambient  fluids  for  mass  conservation 
considerations. 

(5)  For  asymptotic  calculations,  the  entrained  layer  is  thin,  i.  e.,  c  -  R 
or  C  «  1- 

(6)  »  n,  so  Stokes  law  (Batchelor,  1967)  may  be  written  as, 

^  goR^  (1) 


where; 

a  is  the  coefficient  of  thermal  expansion 
g  is  the  acceleration  of  gravity. 
v„  is  the  ambient  kinematic  viscosity 
AT  is  T  - 

(7)  The  portion  of  mass  input  from  the  heater  that  goes  into  the  lengthening 
conduit  is  negligible  (A^U  «  F). 

(8)  The  average  speed  of  fluid  in  the  conduit  (U^)  is  greater  than  U. 

3.4  Thermal  Power  Balance 

Applying  the  principle  of  conservation  of  energy  to  the  stippled  volume  in 
figure  7  yields  the  following  general  relationship  for  thermal  energy  stored  in  the 
system 


where; 
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symbolizes  rate  of  change  per  unit  time 
Oin  is  the  rate  of  heat  transferred  into  the  system 
Oout  is  of  heat  transferred  out  of  the  system 

Q  is  the  thermal  energy  stored  in  the  system. 

Now,  the  heat  input  to  the  system  comes  as  a  result  of  the  warm  fluid  entering 
through  the  circular  region  marked  by  A  (from  the  source)  and  through  the  annular 
surface  region  marked  by  B  (from  entrainment).  Symbolically, 

Qb 

Heat  leaves  the  system  via  conduction  through  the  surfaces  of  the 
hemispheres  marked  by  E  and  D.  So, 

Qo.-Qe-  Qd 

Note  hemisphere  E  has  radius  R  while  hemisphere  D  has  radius  c  (=R  +  C).  The  heat 
transferred  through  E  is  greater  than  that  transferred  through  the  surface  of  the 
hemisphere  marked  by  X,  with  the  difference  being  the  amount  entrained  through 
the  annular  surface  marked  by  B,  i.  e., 

<?,-<?£-  <?r 

Now,  substituting  (3),  (4),  (5),  into  (2)  gives 


—  "  Qa~  Qx  ~Qd 

dt 


(6) 
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Assumption  (1)  allows  us  to  write 

Qa  * 


(7) 


where:  Cp  is  the  specific  heat  at  constant  pressure. 

AT,=T,  -  (constant) 

To  find  an  expression  for  the  heat  lost  from  the  system  at  the  surface  of  the 
sphere  marked  by  X  and  D  (radius  =  c),  use  assumptions  (2)  and  (3).  Then, 


Qx*  Qd^  1,^ 


.  .  ,e*Ar. 
-  4it*Y(— — ) 


(8) 


where: 

k  is  the  thermal  conductivity 

Y  is  a  dimensionless  parameter,  defined  so  that  is  the  ratio  of  heat 
loss  from  the  sphere  of  radius  c  to  total  heat  loss  from  the  sphere  of 
radius  R.  Note  that  Os  y  ^  i- 

Making  these  substitutions  for  Ojj’  ^nd  Op*  equation  (6)  becomes 


diQ) 

dx 


pcFLT,  -  47iJi:Y(^) 
0 


(9) 


The  thermal  energy,  Q,  stored  in  the  system  is 

Q  -  pc^KAT  (10) 

where  V  is  the  volume  of  the  system. 

Substituting  this  expression  into  (9),  using  assumption  (4),  and  recalling  that  k  « 
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pCpK  gives 


djVLD 

dt 


FUT,  -  4;ticY 


(11) 


where  k  is  thermal  diffusivity. 

Note  that  if  the  heat  loss  term,  *  Qj^  ^  is  zero,  this  equation  reduces  to  the 

energy  balance  used  in  the  Griffiths  and  Campbell  (1990)  model. 


3.5  Mass  Balance 

Using  assumptions  (4)  and  (7),  the  application  of  the  principle  of  conservation 
of  mass  to  the  system  yields 


dt 


F  *  e 


(12) 


where;  e  is  the  rate  of  entrained  fluid  entering  the  system,  and  is  the  area 
of  the  conduit. 

Since  the  speed  of  the  fluid  at  the  equator  of  the  ball  is  U/2  (Happel  and 
Brenner,  1965),  ambient  fluid  is  entrained  at  B  at  the  rate  «  .  nJtl/C.  So,  equation 

(12)  becomes 


dt 


F  ♦  kJJC/C. 


(13) 
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3.6  Simplification 

Wc  can  estimate  the  thickness  of  the  thermal  boundary  layer,  6,  by  (xi)^  ■ 

The  time  that  the  hot  and  cold  fluids  are  in  contact  is  proportional  to  R/U.  So, 


where  k|  is  a  dimensionless  constant. 

In  the  best  model  available,  Griffiths  and  Campbell  (1990)  assumed  that  all  of 
the  heat  added  to  the  system  was  entrained  in  the  ball  of  the  plume.  In  the  present 
notation,  this  is  equivalent  to  setting  C  equal  to  6,  i.  e.,  all  of  the  thermal  boundary 
layer  is  entrained.  The  experimental  study  in  chapter  2  suggested  that  not  all  of  the 
heat  added  to  the  system  was  entrained  in  the  system.  This  hcatloss  is  most  vivid  in 
the  T^b  =  -26.1C  case  (figure  6b).  Heat  loss  will  also,  however,  play  a  role  in  the 
other  two  cases  considered.  In  all  of  these  cases,  only  a  portion  of  the  thermal 
boundary  layer  will  be  swept  into  the  ball,  i.  e.,  C  <  6-  Define  a  ratio  between  C  and 
6  using  a  dimensionless  parameter  n  where 

C  =  t,6  ^^5) 

Note  0sT)slsoT]  =  0  means  no  entrainment  and  *  1  means  no  heat  loss. 

This  new  parameter,  n.  is  related  to  the  location  of  the  rear  stagnation  point 
in  figure  7.  The  streamline  around  the  ball  that  ends  at  this  stagnation  point  is  the 
dividing  line  between  entrained  fluid  and  that  left  behind.  The  position  of  this 
stagnation  point  will  vary  depending  on  the  size  of  U  relative  to  the  speed  of  the 
fluid  coming  up  the  conduit  (U^).  A  reasonable  relationship  that  incorporates  this 
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speed  dependence  is 


Tj-1 


(16) 


where: 

P  is  a  nonnegative  dimensionless  constant 

(=  F/Aj)  is  the  average  speed  of  fluid  in  the  conduit. 

For  a  continually  fed  plume,  U  <  U^,  so  t)  <  1,  as  required.  The  constant  p  provides  a 
measure  of  the  effect  that  the  U/U^  ratio  has  on  n-  As  p  approaches  infinity,  n 
approaches  1  and  we  again  have  the  case  in  which  all  of  the  thermal  layer  is 
entrained  (Griffiths  and  Campbell,  1990).  For  p  =  0,  t)  =  0  and  none  of  the  heated 
ambient  fluid  is  entrained  into  the  ball.  For  our  purposes,  the  appropriate  value  of 
P  will  be  determined  from  experimental  data  or  numerical  computation  data. 

The  previously  defined  dimensionless  quantity,  y.  in  equation  (8)  can  be 
expressed  in  terms  of  t).  Since  y^  is  the  ratio  of  heat  lost  from  the  system  at  the 
surface  of  the  sphere  of  radius  c  to  the  total  heat  lost  through  the  surface  of  the 
sphere  of  radius  R  (figure  8),  we  have 


H 


(17) 


Using  Stokes  law  (1),  the  thermal  boundary  layer  approximations  (14)  and  the 


Figure  8.  Simplified  temperature  profile  through  the  thermal  boundary  layer. 


thermal/fluid  boundary  layer  relationships  (15-17),  and  assumption  (S)  yields: 


mass: 


dV 


-  F  * 


U 


(18) 


energy: 


dt 


FLTj 


I2nx}^vj  u  yWi 


09) 
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Nondimcnsionalizc  equations  (18)  and  (19)  by  setting 


,  R  ^  t  a  AT  U  V 

r»— ,  t— ,  0« - ,  — ,  V- — 

R.  I.  Ar.  t/.  K. 


where. 


9v  F  ,  K 

F.-  - - = -  .  K.-4/3nFj.  r.-— 

SngaAT,  "  "  F 


R*  is  the  asymptotic  radius  of  the  ball  for  the  case  of  no  entrainment  or  no  heat  loss 
(Whitehead  and  Luther,  1975).  In  that  case,  the  radius  of  the  ball  will  grow  no  larger 
than  R*.  V*  is  the  volume  of  a  sphere  of  radius  R«  and  t,  is  the  time  required  to 
fill  this  volume  at  a  flow  rate  of  F.  Using  assumption  (7),  the  resulting 
dimensionless  equations  are 


(a)  (b) 


W(1  -  «<») 


d(yQ)  ^ 

(«)  (J)  (g) 


where, 


(  X*  j  V.K* 


27;*/V  9  J  "  k* 
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Note  S  has  been  defined  to  be  analogous  to  the  standard  Rayleigh  number.  The 
Stokes  equation  (1)  in  dimensionless  form  is 


u 


(25) 


where. 


1  j 

0.011 

1  9 

1 

(26) 


At  small  times,  before  ball  liftoff,  the  entrainment  term  (c)  in  equation  (22) 
and  the  heat  loss  term  (g)  of  equation  (23)  can  be  neglected.  With  these 
simplifications,  these  equations  yield  a  rise  law  identical  to  that  of  Whitehead  and 
Luther  (1975): 


e  -  1 

r  -  T*'’  (27) 


These  relationships  hold  for  t  <  tj,  where 

T,  -  (28) 

For  T  >  Xj,  the  entrainment  term  (c)  becomes  important  in  the  mass  equation  (22).  In 
this  case,  the  rise  law  is  identical  to  Griffiths  and  Campbell  (1990): 

r  -  x>^ 

0  -  x-^* 


(29) 
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This  result  will  hold  as  long  as  the  heat  loss  term  (g)  is  negligible,  i.e.,  for  Tj  <  t  < 
Tj,  where 


'^3/p  -  •  1\ 

I 


When  T  >  Tj,  heat  loss  (g)  balances  input  (f)  in  the  energy  equation  (23)  and 
entrainment  (c)  balances  the  rate  of  storage  (a)  in  the  mass  equation  (22),  giving  a 
new  result 


Mg 

r  - 

5«4D 

e  - 

1 

U  -  T^*PP. 


(31) 


3.8  Comparison  with  Experimental  Result 
The  asymptotic  model  developed  in  this  chapter  modifies  the  Griffiths  and 
Campbell  (1990)  model  by  adding  the  effects  of  heat  lost  from  the  ball.  This 
addition  results  in  a  smaller  asymptotic  rise  speed  (equation  31)  for  any  finite  p, 
compared  to  the  best  previous  model  (Griffiths  and  Campbell,  1990).  Since  the  speed 
of  the  experimental  results  was  roughly  constant,  and  the  Griffiths  and  Campbell 
(1990)  rise  law  predicted  slightly  increasing  speeds,  any  finite  choice  of  P  will  bring 
the  new  model  closer  to  the  experimental  results. 

This  improvement  may  be  measured  quantitatively  by  computing  the 
percentage  of  the  difference  between  the  Whitehead  and  Luther  (1975)  model  and 
the  optimal  exponential  model.  This  same  comparison  was  used  to  evaluate  the 
performance  of  the  Olson  and  Singer  (1985)  and  the  Griffiths  and  Campbell  (1990) 


42 


models  in  chapter  2  where  the  best  percentage  improvement  (Griffiths  and  Campbell, 
1990)  was  81.6%  for  the  =  25C  data  and  89.5%  for  the  T^^b  =  case.  The 
new  asymptotic  model  allows  p  to  be  chosen  to  match  the  optimal  exponential  model, 
so  for  p  =  5.76  (for  T^b  =  25C)  and  p  =  1.33  (for  T^b  =  0-lC),  the  percentage 
improvement  is  100%.  The  rise  laws  for  these  cases  are  plotted  in  figures  9  and  10. 
(The  corresponding  volume  versus  time  plots  are  given  in  figures  11  and  12.)  These 
choices  of  P  indicate  a  measurable  improvement,  however,  they  also  suggest  that  the 
addition  of  heat  loss  to  the  model  is  not  as  significant  a  contribution  as  the  addition 
of  entrainment.  The  addition  of  entrainment  alone  (Griffiths  and  Campbell,  1990) 
accounts  for  the  majority  of  the  difference  from  the  experimental  results.  For  the 
Tjujb  =  -26. 1C  case,  p  must  be  chosen  as  -1.866  in  order  to  match  the  optimal 
exponential  model.  This  choice  of  a  negative  p  violates  our  assumption  that  p  >  0 
and  makes  the  new  model  invalid  for  this  case. 

The  determination  of  p  may  be  based  on  comparison  with  experimental  or 
numerical  data,  or  on  some  knowledge  of  the  relationship  between  the  thermal  and 
fluid  boundary  layers  for  a  given  set  of  parameters.  Any  positive,  finite  choice  of  P 
will  improve  the  model.  For  example,  p  =  1  results  in  a  95%  improvement  for  the 
^amb  ®  improvement  in  the  T^b  =  ^.IC  case.  This  general  improvement 

in  these  two  cases  for  all  positive,  finite  p  confirms  that  the  heat  loss  effect  is  an 
important  refinement  of  the  model  when  the  plume  has  the  balloon-on-a-string 
structure. 


3.9  Limits  of  Applicabilitv 

In  deriving  the  new  rise  law  (31),  we  made  several  simplifying  assumptions. 
The  new  rise  law  will  become  invalid  when  these  assumptions  no  longer  reflect  the 
physical  situation.  Two  of  these  assumptions,  one  implicit  and  one  explicit,  appear 
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Figure  11.  Volume  comparison  of  new  rise  law  and  experiment, 
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Figure  12.  Volume  comparison  of  new  rise  law  and  experiment,  T  u-O.lC 


47 


to  present  the  most  severe  limitations. 

3.9.1  Negligible  Conduit  Heat  Loss 

By  neglecting  the  heat  lost  in  the  conduit  compared  to  the  heat  lost  in  the 
ball,  we  assumed  that  the  temperature  of  the  fluid  entering  the  system  through  the 
circular  surface  A  in  figure  7  was  the  same  as  the  source  temperature  of  the  heater. 
As  the  conduit  grows  in  length,  this  assumption  is  no  longer  valid. 

To  estimate  the  height  at  which  the  conduit  length  may  become  important, 
consider  the  approximate  thermal  diffusive  thickness. 


8  = 


For  the  conduit,  the  time  scale  for  this  diffusion  is  U^/z,  where  z  is  the  length  of  the 
conduit.  So,  the  thermal  diffusive  thickness  becomes 


Now,  as  a  first  approximation,  the  conduit  heat  loss  will  become  important  when  this 
diffusive  thickness  equals  the  conduit  radius,  r^.  Recalling  that 


We  have  that 


• 


_F 

Kn 


where  z^  is  the  height  at  which  conduit  heat  loss  is  important. 

For  the  =  25C  case,  z^  =  25.2  cm  and  for  the  =  O.IC  case,  z^  =  18.8 
cm.  This  scaling  argument,  of  course,  only  gives  a  rough  estimate  of  z„,  but  these 
calculations  indicate  that  heat  lost  through  the  conduit  may  be  significant  even  for 
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the  heights  observed  in  the  experimental  studies. 

3.9.2  Balloon-On-A-Strine  Shape 

In  describing  the  system  used  for  the  analysis,  we  implicitly  assumed  that  the 
plume  structure  could  be  modelled  as  a  sphere  trailing  a  narrow  conduit.  The  = 
•26. 1C  case  did  not  exhibit  this  structure.  As  shown  in  section  3.8,  the  new  rise  law 
did  not  apply  to  this  case  since  the  appropriate  choice  of  p  was  negative,  violating 
the  assumption  of  nonnegative  p. 

To  help  understand  some  of  the  factors  that  determine  when  the  balloon-on-a- 
string  shape  will  occur,  consider  the  theoretical  analog  of  a  sphere  in  an  infinite 
tank  of  syrup.  The  sphere  has  radius  a  and  maintains  a  vacuum  internally  (figure 
13).  The  sphere  is  also  able  to  immediately  heat  the  syrup  and  reduce  its  viscosity  to 
zero. 

At  the  surface  of  the  sphere,  the  normal  stress  must  be  continuous,  so  o„  =  -P 
+  pUj  (where  o„  is  the  normal  stress)  is  continuous.  Inside  the  sphere,  however,  P  =  p 
=  0,  so  we  must  have 

AP  =  pa^. 

This  idealized  situation  is  analogous  to  the  heater  in  the  tank.  For  the 
experiment,  then,  this  AP  roughly  represents  the  pressure  difference  needed  to  pull 
cold  syrup  onto  the  heater  surface.  If  this  pressure  difference  is  not  available,  the 
heater  will  draw  preheated  syrup  rather  than  ambient  syrup. 

The  total  pressure  head  available  in  the  plume  is  roughly  ApH,  where  H  is  the 
total  height  of  the  plume.  This  total  must  supply  the  head  to  send  fluid  up  the 
conduit  and  to  draw  ambient  fluid  onto  the  heater.  While  the  ball  is  forming  on  the 
heater,  the  pressure  head  is  increased  by  increasing  the  size  of  the  ball.  When  the 
ball  radius  is  large  enough  to  generate  the  pressure  difference  needed  to  draw 
ambient  fluid,  the  ball  lifts  away  from  the  heater  and  trails  the  narrow  conduit.  In 
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the  =  -26. 1C  ease,  ^lu,  is  so  large  that  the  ball  is  never  able  to  generate  enough 
head  to  pull  in  new  ambient  fluid,  so  it  recycles  previously  heated  fluid. 

Using  the  idealized  spherical  geometry  analog,  we  can  estimate  the  height 
needed  for  liftoff.  Since 


AP  -  ApgH^  = 


and 


F  * 


we  can  solve  for  the  height  needed  for  liftoff. 


In  A  pgr ^ 


where  rj,  is  the  heater  radius 

Hq  is  the  height  needed  for  sphere  liftoff. 

Using  this  relationship  with  the  experimental  parameters,  we  have 

^.lisc  "  -02  cm 


For  the  experiment,  liftoff  for  the  ■  25C  case  occurred  at  approximately  3.5 
cm.  For  the  =  0.1  C  case,  liftoff  occurred  at  approximately  13  cm.  No  liftoff 
occurred  for  the  »  -26. 1C  case. 


CHAPTER  4 


NUMERICAL  MODELLING 

4.1  Purpose 

An  accurate  numerical  model  of  starting  plume  development  provides  several 
significant  advantages  over  the  experimental  investigation.  First,  a  numerical  model 
allows  the  investigation  of  parameter  spaces  that  were  not  or  could  not  be  explored 
in  laboratory  experiments.  Such  models  also  allow  more  detailed  observations  than 
laboratory  experiments.  For  example,  while  velocity  values  may  only  be  observed  at 
points  where  instrumentation  is  available  for  the  experiment,  the  velocity  is  known 
at  each  computational  node  in  the  numerical  model.  There  are,  in  general,  many 
computational  nodes  throughout  the  domain  of  interest,  and  they  are  more  dense  in 
regions  of  rapid  flow  change.  These  nodes,  therefore,  serve  as  a  network  of  ’sensors’ 
that  are  distributed  much  better  than  usually  allowed  in  the  laboratory  environment. 

Since  a  good  numerical  model  will  provide  the  means  to  run  more  detailed 
’experiments’  in  a  wider  parameter  range  than  the  laboratory  case,  verification  of 
the  numerical  model  is  essential.  The  experimental  study  of  chapter  2  provides  an 
excellent  opportunity  to  perform  such  verification.  The  goal  of  the  numerical  model 
in  this  study  is  to  accurately  predict  the  observed  behavior  of  the  experiments, 
principally  the  height  of  the  plume  as  a  function  of  time.  In  subsequent  studies  it 
may  be  used  to  predict  creeping  thermal  plume  behavior  in  other  cases  with  some 
confidence. 
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4.2  Governing  Equations 

The  physical  problem  is  well-described  by  the  continuity  equation  for 
incompressible  fluids,  the  Navier-Stokes  equations  and  the  energy  equation,  with  the 
Boussinesq  approximation.  These  equations  arc  given: 


du. 

^-0  (32) 

at, 


du.  dU;  at.. 


'cbc, 


ar. 

; 


(33) 


with. 


du.  du, 


dT  dr 
pc, —  ♦pc.u, —  ♦ — i=0 
'dr  ^  ^  ^dxj  dxj 


(34) 


(35) 


with, 


(36) 


and  p  varies  with  temperature  as  indicated  in  figure  2  of  chapter  2,  i.  c., 
logio  p  »  20.606  -  14.991(1000/T)  +  3.0762(1 000/T)* 
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Where  n  has  units  mPa-sec  and  T  is  in  Kelvin  units.  In  equations  (32)  through  (36)  t 
is  the  time,  Uj  is  the  velocity  component  in  the  Xj  coordinate  direction,  P  is  the  non¬ 
hydrostatic  pressure,  T  is  the  temperature,  p  is  the  density,  tjj  is  the  stress  tensor,  Qj 
is  the  heat  flux  vector,  p  is  the  dynamic  viscosity,  Cp  is  the  heat  capacity,  k  the 
thermal  conductivity  and  a  is  the  coefficient  of  thermal  expansion.  is  a 
reference  temperature  for  which  buoyant  forces  are  zero,  is  the  identity  tensor, 
and  gj  is  the  gravitational  acceleration. 

4.3  Boundary  Conditions 

The  computational  domain  and  boundary  conditions  have  been  chosen  to 
model  the  physical  situation  in  the  experimental  apparatus.  The  axisymmetry  of  the 
plume  suggests  making  the  domain  axisymmetric  with  the  radial  coordinate 
measured  out  from  the  center  of  the  heater.  Boundary  conditions  include  the  no-sUp 
condition  along  the  top,  side,  and  bottom  walls.  (Since  a  relatively  thick  layer  of 
dried  syrup  formed  on  the  surface  of  the  fluid,  the  boundary  conditions  are  set 
assuming  that  the  vessel  is  completely  enclosed.)  Along  the  axis  above  the  heater, 
normal  velocity  and  heat  flux  are  zero  due  to  symmetry  considerations.  Along  the 
bottom,  perfect  insulation  is  assumed,  excluding  the  heater.  The  temperature  is 
assumed  to  be  ambient  on  the  top  and  side.  The  heater  input  is  specified  as  a 
boundary  condition  on  the  bottom,  corresponding  to  the  physical  size  of  the  heater 
in  the  experiment.  The  value  of  heat  flux  was  calculated  from  the  resistance  of  the 
heater  and  the  applied  voltage  across  the  heater.  Figure  14  presents  a  graphical 
interpretation  of  the  domain  and  boundary  conditions. 

4.4  Computational  Method 


The  most  common  tools  for  solving  incompressible  flow  problems,  such  as  that 
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Figure  14.  Computational  domain  and  boundary  conditions. 
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described  in  the  previous  section,  are  the  finite-difference  method  (FDM)  and  the 
finite-element  method  (FEM).  The  technique  chosen  to  construct  the  numerical 
model  for  the  plumes  was  the  FEM.  The  FEM  has  been  the  method  of  choice  in  the 
computational  fluid  dynamics  literature  for  the  last  few  years.  The  major  reasons 
that  FEM  was  chosen  for  this  problem  were  the  flexibility  it  provides  in  choosing 
the  computational  grid  and  the  availability  of  an  existing  code  that  was  well-suited 
to  the  plume  problem.  The  FEM  allows  a  variable  mesh  to  be  used  when  discretizing 
the  domain;  the  FDM  requires  that  the  nodal  spacing  be  equal  throughout  the 
domain.  This  selection  flexibility  provides  the  opportunity  to  make  the  mesh  dense 
in  areas  where  large  gradients  occur.  For  this  specific  problem,  the  finer  mesh  may 
be  used  to  resolve  the  large  temperature  gradients  on  the  boundary  of  the  plume. 

The  finite-element  code  used  to  produce  the  numerical  model  is  named  NACHOSII. 
NACHOSII,  written  by  D.  Gartling  (1986)  at  Sandia  National  Laboratories,  is 
designed  for  two-dimensional  planar  or  axisymmetric  analysis  of  viscous, 
incompressible  flows,  including  the  effects  of  buoyancy  and  heat  transfer.  A  more 
complete  description  of  NACHOSII  will  be  given  in  section  4.9.1. 

The  brief  description  of  the  FEM  presented  here  is  practically  and 
computationally  oriented.  A  large  body  of  literature  describes  the  more  theoretical 
aspects  of  the  technique.  Two  excellent  references  that  address  the  FEM  from 
theory  and  application  to  fluid  flow  problems  are  Gunzburger  (1989)  and  Cuvelier, 
et.  al.  (1986).  The  following  discussion,  based  largely  on  Gartling  (1986),  describes 
the  FEM  in  general  terms;  specific  application  to  the  plume  problem  is  given  in 
section  4.10. 

The  FEM  begins  by  subdividing  the  domain  of  interest  into  a  number  of 
geometrically  simple  regions  called  finite  elements.  Within  each  element,  the  desired 
unknowns  (uj,  P,  T)  are  located  at  specified  points  called  nodes.  These  unknowns  are 
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interpolated  by  some  suitable  interpolation  function  at  the  nodes.  The  proper  choice 
of  the  elements  and  interpolation  functions  is  crucial  to  the  success  of  the  FEM. 
Requirements  for  this  choice  will  be  considered  in  a  later  section.  Using  the 
notation  of  Gartling  (1986),  the  unknowns,  Qj,  P,  and  T  at  a  general  point  Xj  within 
an  element  may  be  written  as  a  vector  product  of  the  interpolation  functions  <t, 
and  and  the  corresponding  nodal  values: 

«/V)  =  *\x)uj(t) 

7Tx„0  =  e^(x.)7t0. 

The  superscript  T  denotes  transpose. 

Substituting  these  expressions  into  the  governing  equations  (32)  -  (36)  yields 
the  following  set  of  equations. 

Momentum: 

/^($.f  .e,u..F.7)  = 

Continuity. 

=  R, 

Energy: 

=  Rj. 

These  equations  will,  of  course,  no  longer  be  an  exact  representation  of  the  physical 
interactions  described  by  (32)  -  (36).  The  amount  that  they  differ  due  to  the 
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interpolation  approximation  is  represented  by  the  residuals,  or  errors,  denoted  by 
Itui,  ^T-  FEM  seeks  to  minimize  these  residuals  in  some  weighted  sense. 

The  particular  FEM  used  to  model  the  plumes  is  a  Galcrkin  method.  In  this 
technique,  the  residuals  are  minimized  by  making  them  orthogonal  to  the 
interpolations  functions  over  each  element.  This  condition  may  be  expressed  as 

=  0 

=  0  (37) 

=  (T/p  =  0 

(e,«j)  =  (e/j)  =  0 


where  <.,.>  is  defined  as 

id, Si  =  f^d-bdQ 


with  Q  being  the  area  of  an  element. 

The  Galerkin  method  reduces  the  governing  partial  differential  equations 
(32)  -  (36)  to  a  nonlinear  system  of  ordinary  differential  equations  in  time.  The 
unknowns  of  this  system  are  (ilj,  F,  T)  for  each  element.  The  coefficients  of  the 
system  are  derived  from  the  integrals  represented  by  the  inner  products  in  equation 
(37).  The  details  of  this  manipulation  are  described  in  the  appendix.  The  results  can 
be  expressed  in  matrix  form  as  follows 
momentum: 


(38) 
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I 

I  continuity: 


-(fV=Q 


energy: 


(39) 


N—*IKU)f*IXf)f=Cif,U)  (40) 

dt 


where 


(41) 


Each  of  these  element  matrices  and  vectors  has  a  physical  meaning  associated 
with  the  terms  in  the  governing  equations.  The  M  and  N  matrices  Teprescnt  the  mass 
and  heat  capacity  terms  in  the  momentum  and  energy  equations,  respectively.  C  and 
D  are  the  result  of  the  advection  terms  in  the  momentum  and  energy  equations.  K 
and  L  represent  the  diffusion  terms.  The  Q  matrix  is  the  gradient  operator  while 
is  the  divergence  operator.  The  B  matrix  represents  the  buoyancy  term  from  the 
Boussinesq  approximation.  E  and  <j  allow  for  forcing  functions  in  the  system;  for 
this  problem,  is  zero.  These  physical  interpretations  are  very  helpful  in  analyzing 
the  method’s  performance  for  the  given  plume  application.  The  exact  integral 
expressions  for  these  matrices  are  given  in  the  appendix.  Once  the  matrices  arc 
developed  for  one  element,  the  global  matrix  system  may  be  assembled  from  the 
element  matrices  by  addition,  i.  e.. 
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^glM  ~  ]C 


<•1 


etc. 


where  nelm  is  the  total  number  of  elements  in  the  grid.  Note  that  the  global  system 
has  the  identical  form  of  (38)  -  (41),  where  the  matrices  represent  global  matrices. 
Equations  (38)  -  (41)  may  then  be  combined  into  a  single  matrix  equation 


1 

{dO\ 

/MOO) 
0  0  0 

dt 

dP 

C(U)*K(0.T)  -Q  B(J)  ' 

-Q*  0  0 

ft/] 

P 

0 

,0  0  N, 

dt 

df 

dtj 

0  Q  IKU)  *  Un 

(42) 


or, 


M—  *  K{U,f)V  =  FiU.f) 
dt 


(43) 


where 


There  are  two  common  methods  for  solving  the  matrix  system  (42):  the  mixed  or 
integrated  FEM  and  the  penalty  FEM. 

4.4.1  Integrated  Finite  Element  Method 

The  integrated  FEM  is  a  straightforward  extension  of  the  process  described 
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above.  The  system  of  ordinary  differential  equations  in  time  is  solved  for  velocity, 
pressure  and  temperature  at  each  time  step 
4.4.2  Penalty  Finite  Element  Method 

In  contrast  to  the  integrated  FEM,  the  penalty  FEM  eliminates  pressure  as  an 
unknown  in  the  matrix  system.  Since  a  large  part  of  the  cost  of  any  FEM  is  the 
solution  of  the  matrix  system,  the  penalty  method  generally  has  reduced 
computational  and  storage  costs  when  compared  to  the  integrated  FEM.  The 
theoretical  background  for  the  penalty  method  comes  from  the  application  of 
constrained  optimization  theory  to  a  linearized  form  of  (32)  and  (33).  Let  us  briefly 
consider  this  form  in  order  to  develop  the  equations  used  in  the  penalty  method. 
When  the  advection  and  buoyancy  terms  in  (32)  are  dropped,  the  resulting  linear 
equations  are  commonly  called  the  Stokes  equations 


dUf  dx^ 
dt  ~  dxj 

du. 

— ^  =  0 
dx, 


0 


This  problem  may  be  recast  as  a  minimization  problem  where  the  momentum 
equation  is  constrained  by  the  continuity  equation.  A  standard  method  of  enforcing 
such  a  constraint  is  the  method  of  Lagrange  multipliers  with  penalization  (Girault 
and  Raviart,  1986).  This  approach  explains  the  use  of  the  term  penalty  for  this 
method. 

Returning  now  to  the  nonlinear  problem  (32)  -  (36),  recall  that  it  has  no 
known  equivalent  minimization  form,  but  the  penalty  method  may  be  applied  to  this 
problem  by  using  a  slightly  different  approach.  When  the  minimization  form  of  the 
Stokes  problem  is  solved  by  the  penalty  method,  the  corresponding  Euler-Lagrange 
equations  are  identical  to  the  Stokes  equations,  except  that  the  continuity  equation  is 
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relaxed,  i.  e.. 


dx. 


-€P. 


The  parameter  e  is  the  penalty  parameter.  Note  that  the  continuity  equation  is 
recovered  as  e  goes  to  zero.  One  can  think  of  this  as  a  pertubation  to  the  continuity 
equation  that  introduces  some  small  compressibility  into  the  system. 

Applying  this  perturbation  to  the  continuity  equation  (32),  approximating  the 
velocity  and  pressure  with  the  interpolation  functions,  and  using  the  Galerkin 
procedure  gives 


Or,  in  matrix  form 


Q^U  =  -€M/ 


where 


Ja 


Solving  the  matrix  equation  for  P  gives 

€ 


This  equation  allows  pressure  to  be  eliminated  from  the  global  matrix  system.  The 
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global  system  for  the  penalty  method,  therefore,  may  be  written  as  a  matrix  equation 


[m  o' 

loiv, 


dt 

df 


\dt) 


c(U)^K(u,T)^Kj.  am  Yl/]  _  f  1 

0  DiU)  *  L(7)AfJ  " 


where 

K,  - 

€ 


After  the  system  is  solved  for  the  velocity,  the  pressure  may  be  recovered  using  the 
matrix  form  of  the  perturbed  continuity  equation. 

Note  that  to  allow  for  easy  element-level  construction  and  subsequent  global 
matrix  assembly  of  the  Kp  matrix,  it  is  necessary  that  the  terms  of  Mp**  be  easily 
constructed,  i.  e.,  Mp  must  be  invertible  at  the  element  level.  Since  Mp  is  defined  in 
terms  of  the  pressure  interpolation  functions,  these  functions  may  only  be  defined 
within  an  element,  i.  e.,  they  can  be  discontinuous  between  elements.  For  this  reason, 
the  pressure  interpolation  functions  for  the  penalty  method  will  always  be 
discontinuous  between  elements. 


4.4.3  Iterated  Penalty  Method 

Gunzburger  (1989)  describes  an  iterated  penalty  approach  that  enforces  the 
incompressibility  condition  using  an  iterative  application  of  the  perturbation  of  the 
continuity  equation.  The  algorithm  is  given  as 


1.  Given  solve  for  fj",  T" 
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+  c(y")i/"  -  =  6 

at 

N—  *  DiU*)!*  +  L(f^f  =  6 
dt 


2.  Update 


p"  .  P*-*  +  J-A/ 

e  ' 

(This  step  is  based  on  a  perturbation  to  the  continuity  equation.) 

3.  If  discrete  divergence  is  too  large  repeat.  (The  traditional  penalty  method 
occurs  in  the  method  when  Pq  =  0  and  steps  2  and  3  are  performed  only  once.) 

Results  of  numerical  experiments  with  this  method  were  generally  poorer  than  those 
found  with  the  Uzawa  algorithm  described  later  (4.9.4).  It  is  included  here  as  a 
reference  to  compare  with  the  Uzawa  algorithm. 

4.4.4  Time  Integration 

Both  the  integrated  and  penalty  FEM  provide  a  spatial  discretization  that 
results  in  a  system  of  first-order  ordinary  differential  equations  in  time.  The 
solution  of  the  system  is  found  by  finite  differencing  in  time;  see  Strikwerda  (1989) 
and  Sod  (1984).  In  the  NACHOSII  code,  a  predictor-corrector  method  is  employed 
that  uses  an  explicit  time-differencing  for  the  predictor  and  an  implicit  method  for 
the  corrector.  The  implicit  differencing  is  unconditionally  stable  for  all  choices  of 
time  steps  and  thus  avoids  severe  restrictions  in  the  choice  of  the  time  step.  A  choice 
of  first-order  or  second-order  methods  is  available.  The  first-order  method  uses  a 
forward  Euler  predictor  and  backward  Euler  corrector.  The  second-order  method 
uses  Adams-Bashforth  predictor  and  trapezoid  corrector.  The  second  order  methods 
have  the  advantage  of  increased  accuracy  while  the  first  order  Euler  methods  have 


64 


damping  characteristics  that  are  often  helpful  in  smoothing  the  starting  conditions 
in  the  initial  stages  of  a  problem.  Automatic  time-step  control  is  also  provided  with 
both  methods. 


4.5  Element  Choice 

This  choice  of  elements  for  a  finite-element  model  is  a  matter  of  determining 
the  appropriate  geometry  (e.  g.,  triangular  or  quadrilateral)  and  the  order  of  the 
interpolation  functions  defined  on  the  element.  The  regular,  quadrilateral  shape  of 
the  computational  domain  suggested  the  use  of  quadrilateral  elements.  NACHOSII 
provides  an  element  library  that  includes  an  eight-node  quadrilateral  and  a  nine- 
node  quadrilateral  (figure  15).  These  elements  will  be  considered  specifically  in  the 
folljwing  discussion. 

Furthermore,  the  choice  of  the  order  of  the  interpolation  functions  is  an 
important  consideration.  Cuvelier,  et.  al.  (1986)  describe  necessary  conditions  for 
these  functions.  When  Green’s  formula  is  applied  to  the  Galerkin  form  of  the 
gov  ;rning  equations  (equation  (37)  and  the  appendix),  first-order  derivatives  of  the 
ve)  'City  interpolation  functions  occur,  but  the  pressure  interpolation  functions  arc 
no;  differentiated.  This  means  that  the  velocity  functions  must  be  continuously 
dif  ferentiable  in  each  element  and  continuous  in  the  whole  domain.  The  pressure 
functions,  however,  need  oniy  be  continuous  in  each  element;  they  may  be 
discontinuous  between  elements.  Cuvelier,  ct.  al.  (1986)  also  indicate  that,  in  order  to 
prevent  the  velocity  approximation  from  being  totally  specified  by  the  continuity 
equation,  one  generally  requires  that  the  pressure  interpolation  functions  be  at  least 
one  order  less  than  the  velocity  interpolation  functions. 

NACHOSII  gives  several  interpolation  choices  that  arc  consistent  with  these 


criteria.  In  each  choice,  the  same  interpolation  function  is  used  for  temperature  and 


Figure  15.  Eight  and  nine  node  quadrilateral  elements. 

velocity  components.  For  velocity,  biquadratic,  continuous  interpolation  functions 
may  be  used  for  the  eight-node  quadrilateral  (QUADS)  or  the  nine-node 
quadrilateral  (QUAD9).  For  the  pressure  approximation,  one  may  choose  continuous, 
bilinear  functions  or  discontinuous  linear  functions. 

The  element  and  interpolation  choices  arc  crucial  to  the  success  of  the 
Galerkin  FEM  since  they  determine  which  discrete  finite-clement  subspacc  is  used  to 
approximate  the  solution.  This  choice  is  the  topic  of  much  current  research.  The 
following  discussion  gives  some  basic  theoretical  guidlines.  Gunzburger  (89)  gives  a 
description  of  which  spaces  are  allowed  for  approximating  equations  (32)  -  (36). 
Choosing  velocity,  pressure,  and  temperature  interpolations  arbitrarily  will  often 
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lead  to  poor  results.  In  incompressible  flow  calculations,  the  most  difficult  condition 
to  satisfy  is  the  divergence-free  condition  (32).  Gunzburger  (1989)  also  discusses  a 
key  theoretical  condition  that  relates  to  the  satisfaction  of  the  continuity  equation 
(32).  It  is  known  as  the  Ladyzhenskaya-Babuska-Bressi  (LBB),  the  inf-sup,  or  the 
div-stability  condition.  Girault  and  Raviart  (1986)  give  a  complete  discussion  of  this 
condition,  along  with  various  proofs  for  several  elements.  There  are  several 
equivalent  statements  of  the  div-stability  condition;  proof  of  this  condition  is  not 
trivial  and  will  not  be  attempted  here.  One  important  product  of  this  condition  is 
that  it  allows  us  to  predict  the  convergence  rate  of  the  FEM: 

\u  -  u‘‘\^  *  [P  -  P%  ^  \T  -  T‘‘\^  = 

where 


u  =  iu^^) 
(p^q)  = 

l9lo  =  (9.9)'° 


l9l.  = 


and  the  superscript  d  denotes  the  computed  solution.  The  parameters  k  and  1  are  the 
orders  of  the  velocity  (or  temperature)  and  pressure  interpolations  respectively. 

Note  that  d  is  a  parameter  that  relates  to  the  size  of  the  mesh  spacing.  The  QUAD9 
elements  with  continuous  and  discontinuous  pressure  approximation  satisfy  the  div- 
stability  condition  (Girault  and  Raviart,  1986;  Gunzburger,  1989).  The  QUADS 
elements  with  discontinuous  and  continuous  pressure  approximation  do  not  satisfy 
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the  div-stability  condition  in  the  global  sense.  However,  they  do  satisfy  the 
condition  in  a  local  sense,  and  this  consideration  allows  the  same  rate  of  convergence 
to  be  proven  for  these  elements  as  for  the  continuous  pressure  elements  (Girault  and 
Raviart,  1986).  So,  from  a  rate  of  convergence  point  of  view,  none  of  the 
quadrilateral  elements  has  any  theoretical  advantage. 

Another  consideration  in  choosing  appropriate  elements  is  the  phenomenon 
often  called  mesh  locking.  Roughly  speaking,  in  satisfying  the  continuity  equation 
(32),  some  elements  enforce  too  many  incompressibility  constraints.  When  this 
happens,  the  only  solution  is  the  trivial  (or  locked)  solution.  Elements  that  satisfy 
•the  div-stability  condition  will  not  lock.  The  mathematical  analysis  of  this  condition 
is  beyond  the  scope  of  this  paper.  Hughes  (1986)  presents  a  heuristic  method  to 
determine  an  element’s  propensity  to  lock.  While  this  method  is  not  rigorous,  it 
provides  some  means  of  comparing  the  effectiveness  of  different  elements  to  model 
incompressible  flow.  The  following  paragraph  is  a  summary  of  this  method  of 
constraint  counting. 

On  a  standard  mesh  for  a  two  dimensional  problem,  let  n^^  be  the  total 
number  of  velocity  equations  after  boundary  conditions  have  been  imposed  and  let 
n^  be  the  total  number  of  incompressibility  constraints  as  specified  by  pressure 
nodes.  Define  the  constraint  ratio,  r,  by 


The  idea  behind  the  method  is  that  r  should  mimic  the  behavior  of  the  number  of 
equilibrium  equations  divided  by  the  number  of  incompressibility  conditions  for  the 
governing  system  of  partial  differential  equations  (i.  e.,  the  number  of  space 
dimensions  and  1,  respectively).  In  two  dimensions,  then,  the  ideal  value  for  r  is  2. 
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A  value  less  than  2  would  indicate  a  tendency  to  lock  while  a  value  greater  than  2 
indicates  that  not  enough  incompressibility  constraints  are  present.  A  value  of  r  that 
is  less  than  or  equal  to  1  indicates  more  constraints  than  degrees  of  freedom  and 
severe  locking  is  expected.  For  two-dimensional  problems,  the  following  rules 
summarize  the  method: 

r  >  2  too  few  incompressibility  constraints 
r  =  2  optimal 

r  <  2  too  many  incompressibility  constraints 
r  s  1  locking 

Table  4.  Div-stability  and  mesh  locking  for  elements. 


Element 

r 

LBB  satisfied? 

QUAD  8,  disc  P 

2 

no 

QUAD  8,  cont  P 

6 

no 

QUAD  9,  disc  P 

2  2/3 

yes 

QUAD  9,  cont  P 

8 

yes 

The  div-stability  condition,  roughly  speaking,  insures  that  locking  will  not  occur, 
however,  it  does  not  guarantee  an  optimal  value  of  r.  Note  that  a  value  of  r  greater 
than  two  may  imply  that  the  incompressibility  condition  is  not  adequately  satisfied, 
so  an  element  may  satisfy  the  div-stability  condition  (so  it  won’t  lock),  yet  it  may  not 
enforce  the  incompressibility  condition  adequately.  Table  4  summarizes  the  div- 
stability  and  locking  characteristics  of  the  elements  available  in  NACHOSII.  The  r 
values  indicate  that  none  of  the  elements  being  considered  will  lock;  QUADS 
discontinuous  pressure  has  the  optimal  r  value. 

Pelletier,  et.  al.  (1989)  show  that  a  discontinuous  pressure  approximation  is 
more  effective  in  modelling  the  divergence-free  condition  than  a  continuous 
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pressure  approximation.  Their  argument  is  based  on  element  mass  conservation  and 
the  degrees  of  freedom  required.  This  latter  argument  considers  the  dual  role  of 
pressure.  In  the  integral  form  of  the  continuity  equation  (appendix),  the  pressure 
approximation  must  enforce  the  continuity  equation.  Additionally,  the  pressure 
approximation  must  balance  the  viscous  and  buoyancy  terms  in  the  momentum 
equation.  The  pressure  must  adjust  itself,  then,  to  satisfy  the  continuity  equation 
while  matching  the  needed  balance  in  the  momentum  equation.  There  must  be 
enough  pressure  degrees  of  freedom  to  satisfy  both  of  these  requirements.  One  way 
to  enlarge  the  pressure  space  is  to  move  from  a  continuous  pressure  approximation  to 
a  discontinuous  pressure  approximation. 

4.6  Discrete  Divergence 

To  this  point,  we  have  described  the  FEM  in  general  terms.  We  have 
considered  the  general  procedure,  (he  mixed  and  penalty  methods,  and  the  choice  of 
elements  and  interpolation  functions  without  applying  any  methods  to  the  plume 
problem.  From  this  point  forward,  we  will  focus  on  the  plume  problem  application. 
The  basic  theory  that  we’ve  considered  so  far  provides  some  assistance  in  choosing 
which  method  and  element  is  best  suited  for  this  problem,  but  we  need  additional 
tools  to  determine  which  method  will  give  the  best  results  for  the  plume  application. 
The  experience  of  our  calculations  will  show  that  computing  the  discrete  divergence 
is  a  good  diagnostic  for  making  decisions  for  practical  applications  like  the  plume 
problem. 

As  suggested  by  the  previous  discussion  of  the  div-stability  condition,  one  of 
the  most  difficult  aspects  of  numerically  modelling  incompressible  flow  is  the 
satisfactions  of  the  continuity  equation.  This  difficulty  is  compounded  in  the  plume 
problem  by  the  strong  relationship  between  viscosity  and  temperature.  This 
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relationship  results  in  a  problem  with  strong  coupling  between  the  energy  and 
momentum  equations.  It  also  provides  regions  of  high  viscosity  gradients,  for 
example,  at  the  boundary  of  the  plume. 

Pelletier,  et.  al.  (1989)  proposed  a  simple  diagnostic  tool  for  measuring  the 
effect  of  large  (constant)  viscosity  and  large  viscosity  contrast  on  incompressible 
flows.  They  argued  that  traditional  approaches  often  resulted  in  poor  solutions 
because  of  poor  satisfaction  of  the  continuity  equation  or  the  divergence-free 
condition  (even  when  the  elements  used  satisfied  the  div-stability  condition!),  and 
they  suggested  that  computing  an  approximation  to  the  divergence  of  the  velocity 
field  (i.  e.,  a  ’discrete  divergence’)  is  an  inexpensive  measure  of  the  reliability  of  the 
solution. 

In  order  to  study  the  effects  of  viscosity  variation  on  the  plume  problem, 
NACHOSII  was  modified  to  report  the  discrete  divergence  at  each  time  step. 
Following  Pelletier,  et.  al.  (1989),  the  discrete  divergence  was  computed  as 

Z)i3/P=max^^  j  Q^U ^ 

where  is  the  velocity  vector  on  an  individual  clement,  N.  The  maximum 
is  taken  over  all  elements. 

This  modification  was  easy  to  implement,  and  relatively  inexpensive  to 
compute,  since  for  each  element,  the  matrix  equivalent  of  the  divergence  operator, 
Q^,  was  already  computed.  The  product  Q’^0  provided  the  ’components’  of  the 
discrete  divergence;  the  maximum  of  the  absolute  values  of  these  components  was 
chosen  for  the  discrete  divergence.  FEM  theory  predicts  that  the  discrete  divergence 
should  be  on  the  order  of  machine  zero  (10'^®  for  the  Cray-YMP  used  for  the 
calculations)  for  integrated  FEM  and  on  0(e)  for  penalty  FEM  for  a  solution 
demonstrating  acceptable  mass  conservation. 
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In  practical  computations.the  discrete  divergence  served  as  a  good  indicator 
of  the  reliability  of  the  solution.  In  tests  of  the  mixed  and  penalty  FEM,  the  discrete 
divergence  would  generally  increase  by  one  or  two  orders  of  magnitude  before  a 
problem  occurred  in  the  solution.  This  problem  may  be  a  matrix  ill-condiiioning 
error  or  an  unrealistic  value  of  velocity  or  temperature.  This  behavior  of  the 
discrete  divergence  was  confirmed  for  constant  viscosity  and  variable  viscosity 
applications.  Experience  has  shown  that  the  change  in  the  discrete  divergence  is  a 
conservative  diagnostic,  i.  e.,  an  increase  in  discrete  divergence  does  not  always 
signal  a  problem.  It  does,  however,  provide  an  inexpensive  means  of  testing  a 
solution  and  we  will  use  it  as  an  indicator  of  which  scaling,  element,  and  method  is 
best  suited  for  modelling  the  plume  problem. 

4.7  Scaling 

Using  a  dimensionless  form  of  equations  (32)  -  (36)  has  several  advantages 
over  the  dimensional  form.  One  advantage  is  that  such  scaling  can  provide  an 
indication  of  the  relative  importance  of  terms  in  the  equations.  Dimensionless  forms 
can  also  reduce  the  differences  in  the  magnitudes  of  different  terms  in  the 
equations.  Most  related  numerical  studies  (table  1)  cast  equations  (32)  -  (36)  in 
Prandtl  number-Rayieigh  number  form  by  making  the  following  substitutions: 


f-HL  r'-l 

h  K 


(44) 


T 


(45) 
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where  ATg  is  Theater '  ^amb’  ^  height  of  tank,  is  the  viscosity  of  the  ambient 
fluid.  These  substitutions  put  the  equations  in  the  following  form  (dropping  primes 
for  nondimensional  variables): 


uu. 

—  =0 
3c. 


Pr\  dt 


du-  dr,.. 


at  ^ 

j 


with. 


du.  dUi 

T  =-F6.+u(— i— 

«  *  '^^ar,  dx. 


This  Pr-Ra  scaling  is  commonly  used  for  modelling  Benard-type  problems  for  a  thin 
layer  of  fluid.  In  these  types  of  problems,  one  may  study  the  effects  of  increasing 
buoyancy  by  increasing  the  Rayleigh  number.  One  generally  associates  the  Rayleigh 
number  with  the  temperature  difference  across  the  thin  layer,  or,  alternatively,  with 
the  heat  input  to  the  thin  layer.  When  this  type  of  scaling  is  applied  to  the 
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experiment  described  in  chapter  2,  the  Rayleigh  number  is  very  large  due  to  the 
large  height  of  the  tank.  For  =  25C,  the  Rayleigh  number  is  6.5  X  10®  and  for 
^amb  =  0-lC,  the  Rayleigh  number  is  8.5  X  10\  These  large  values  of  Ra  place  too 
much  weight  on  the  buoyancy  term  in  the  momentum  equation  in  the  Pr-Ra  scaling. 

This  imbalance  potentially  creates  rounding  errors  in  the  solution  of  the  FEM 
matrix  equations  which  result  in  a  discrete  divergence  that  is  too  large.  The  global 
matrix  terms  corresponding  to  the  buoyancy  force  (terms  of  matrix  B  in  equation 
(42))  arc  much  larger  than  those  of  the  pressure  and  viscous  terms  (terms  in  matrices 
Q  and  K,  respectively,  in  equation  (42)).  For  some  cases,  this  size  difference 
measured  as  large  as  10®.  These  large  discrepancies  create  the  possibility  of  round¬ 
off  errors  and  poor  matrix  conditioning.  In  fact,  these  errors  accumulate  quickly, 
preventing  successful  time  integration.  As  a  result,  no  successful  Pr-Ra  scaled  cases 
ran  to  completion  despite  numerous  attempts  using  many  different  grids.  Table  5 
illustrates  the  effect  on  discrete  divergence  of  increasing  Ra  on  an  otherwise 
unchanged  problem.  Note  that  the  values  in  this  table  were  for  a  mixed  method 
calculation,  but  similar  results  hold  for  penalty  methods.  The  increase  in  discrete 
divergence  implies  that  mass  is  not  conserved  in  the  solution  for  larger  values  of  Ra. 
Table  5.  Effect  of  Ra  on  discrete  divergence. 


Ra 

Discrete  Divergence 

1.0  X  10^ 

-10'® 

1.0  X  10® 

-10-* 

1.0  X  10® 

-10-® 

One  alternative  to  correct  this  scaling  imbalance  is  to  consider  a  thinner  layer 
of  fluid.  For  example,  one  may  study  the  plume  formation  and  travel  in  the  bottom 
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fourth  of  the  tank.  A  similar  approach  was  used  by  Olson,  ct.  al.  (1987)  for  studying 
plume  formation.  One  goal,  however,  of  this  numerical  study  was  to  test  the  rise 
laws  found  in  chapter  3  for  the  full  plume  flight  path,  and  to  consider  the 
effectiveness  of  the  rise  laws  in  other  geometries  in  later  studies.  Restricting  the 
depth  of  fluid  would  severely  hamper  the  utility  of  the  numerical  model  for  studies 
of  developed  plumes. 

The  approach  used,  instead,  for  the  numerical  model  in  this  study  was  to 
rescale  equations  (32)  -  (36)  to  balance  the  appropriate  terms  in  the  equations  based 
on  the  important  physical  interactions  in  the  problem.  During  the  ball  formation 
and  rise,  the  key  balance  of  forces  in  the  momentum  equation  (33)  occurs  between 
the  buoyancy,  pressure,  and  viscous  forces.  The  inertia  terms  will  have  a  relatively 
small  effect.  With  this  balance  in  mind,  nondimensionalize  the  governing  equations 
with  the  following  substitutions; 


It'  -  —  •  t'  -  —  ■  '  -  — 

"  °  L  "  L 

-il  ;  =  JL  ;  />/  =  _f_ 
AT,  |i^  pgaLTL 


where 


pgaAT, 


K 


AT. 


/  mmtf 


-  T. 


These  substitutions  put  the  equations  in  the  following  form  (dropping  primes  for 
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nondimensional  variables) 


=0 


(46) 


1  du,  du.  dz,i 


(47) 


dT  ST 
— *u. - «0. 

A  ^Sx^  Sx  ? 


(48) 


with, 


du.  du. 


(49) 


This  dimensionless  form  shows  the  proper  balance  between  pressure,  viscous  forces 
and  buoyant  forces.  Note  that  this  dimensionless  form  may  be  readily  implemented 
on  any  primitive-variable  finite-element  code  (such  as  NACHOSII)  by  setting 
thermal  conductivity  and  acceleration  of  gravity  to  1,  heat  capacity  and  coefficient 
of  thermal  expansion  to  Pr,  and  density  to  1/Pr. 

This  scaling  is  equivalent  to  matrix  rescaling  that  has  been  used  successfully 
in  many  applications,  e.  g.,  Koch  (1985)  used  matrix  scaling  to  reduce  the  condition 
number  of  relevant  matrices  by  as  much  as  five  orders  of  magnitude.  For  a 
simplified  example,  suppose  we  have  the  following  matrix  system: 
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K  -Q  B 

L  . 

U 

F 

-Q^  0  0 

P 

= 

0 

0  0  D 

T 

G 

where  K  and  B  have  elements  of  greatly  different  sizes.  To  rescale  them  to  the  same 
order,  let 


m  =  ixBi 


which  is  the  same  as  the  change  of  variable 


B* 

P* 

r 


kB 

p 

X 

1 

x' 


This  rescaling  results  in  a  system  that  has  elements  whose  sizes  are  comparable  which 
can  be  solved  with  less  round-off  error  than  the  original  system 


AT*  -<?•  B' 
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F 
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o 
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T\ 

G 

The  new  scaling  of  equations  (32)  -(36)  resulted  in  marked  improvement  in 


the  performance  of  the  FEM  code  over  that  with  the  Pr-Ra  scaling  by  forcing  the 
terms  in  the  matrices  representing  pressure,  buoyancy,  and  viscosity  to  be  roughly 
the  same  size  (Table  6).  This  scaling,  therefore,  minimized  the  possibility  of  round- 
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Table  6.  Typical  matrix  norms  for  different  scalings;  mixed  FEM,  discontinuous  P, 
Ra=6.5  X  10®. 


Scaling 

Diffusion 

Matrix(K) 

Buoyancy 

Matrix(B) 

Discrete 

Divergence 

Pr  -  Ra 

-10-* 

-10* 

-10*10 

new  scaling 

-10* 

-10* 

-10'*^ 

off  errors  and  provided  an  effective  means  of  computing  the  solution. 


4.8  Effects  of  Element  Choice 

When  we  discussed  elements  previously,  we  gave  some  criteria  for  choice 
based  on  the  div-stability  condition  and  the  heuristic  argument  of  mesh  locking. 
From  these  criteria,  we  can  conclude  that  all  of  the  available  elements  have  the  same 
theoretical  convergence  properties,  but  that  the  elements  with  discontinuous  pressure 
interpolation  functions  have  better  divergence-constraint  properties. 

The  discrete  divergence  calculations  confirm  these  characteristics  (table  7). 
Note  that  these  results  were  calculated  using  the  mixed  method  (since  penalty 
method  requires  discontinuous  pressure  interpolation)  on  a  sample  grid  for  the  T^j, 
=  25C  conditions.  In  both  QUADS  and  QUAD9  elements,  the  discontinuous  pressure 
approximation  had  a  smaller  discrete  divergence  than  the  continuous  pressure 
approximation.  Of  these  two  elements,  the  one  with  the  optimal  value  of  r  has  the 
smaller  discrete  divergence.  Note  that  this  element  does  not  satisfy  the  div-stability 
condition  in  the  global  sense,  so  while  this  condition  is  related  to  mass  conservation, 
it  is  not  sufficient  to  insure  incompressibility.  These  results,  as  well  as  similar 
results  in  other  ambient  conditions,  suggest  the  best  available  element  to  use  for  the 
plume  calculation  is  the  QUADS  with  linear,  discontinuous  pressure  approximation. 
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Table  7.  Effect  of  element  choice  on  discrete  divergence. 


Element 

r 

Discrete  Divergence 

QUAD  8,  disc  P 

2 

_10-i5 

QUAD  8,  cont  P 

6 

-10-’ 

QUAD  9,  disc  P 

2  2/3 

QUAD  9,  cont  P 

8 

-10*’ 

4.9  Solution  Techniques 

The  choice  of  scaling  and  elements  that  was  discussed  in  earlier  sections  kept 
the  discrete  divergence  to  a  minimum  for  the  sample  computations  performed  for 
comparison.  As  the  solution  progressed  in  time,  however,  the  discrete  divergence 
would  often  grow.  Gradual  growth  was  expected  due  to  the  larger  values  of  speed 
and  accumulation  of  round-off  error,  but  rapid  growth  served  as  a  good  indicator 
that  the  solution  would  encounter  an  ill-conditioned  matrix  error,  or  that  the 
computed  temperature  and  velocity  would  increase  to  unrealistically  large  values. 

Four  different  methods  were  used  in  attempting  to  solve  the  plume  problem. 
The  first  two  were  already  coded  into  the  NACHOSII  system;  the  code  was  modified 
to  accomodate  two  additional  methods  that  were  variations  of  the  first  two  methods: 

1.  The  integrated  or  mixed  method 

2.  The  penalty  method 

3.  The  Uzawa  algorithm. 

4.  The  mixed  method  using  a  multistep  Newton  method 


The  following  paragraphs  will  report  on  the  results  using  these  methods.  For 


each  method,  many  different  tests  involving  various  ambient  conditions  and 
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assumptions  were  performed.  Additionally,  several  meshes  and  time  steps  were  used 
for  each  problem.  Before  presenting  the  results  for  each  of  these  methods,  we  will 
briefly  describe  the  NACHOSII  system. 

4.9.1  The  NACHOSII  System 

The  finite-element  code  used  to  produce  the  numerical  model  is  named 
NACHOSII.  NACHOSII,  written  by  D.  Gartling  (1986)  at  Sandia  National 
Laboratories,  is  designed  for  planar  or  axisymmetric  two-dimensional  analysis  of 
viscous,  incompressible  flows,  including  the  effects  of  buoyancy  and  heat  transfer. 
NACHOSII  represents  a  significant  improvement  over  its  predecessor,  NACHOSI. 
Enhancements  include  a  wider  choice  of  elements,  interpolation  functions,  and  time 
solution  procedures.  This  flexibility  provided  a  combination  of  factors  that 
produced  a  more  effective  numerical  model  than  that  allowed  by  NACHOSI. 

Another  key  difference  between  NACHOSI  and  NACHOSII  is  that  NACHOSII  solves 
the  momentum  and  energy  equations  as  one  system.  NACHOSI  solves  each  equation 
individually  and  then  uses  the  solution  to  solve  the  remaining  equation. 

NACHOSII  provides  a  choice  of  triangular  or  quadrilateral  elements.  The 
velocity  and  temperature  are  approximated  by  quadratic  functions  that  arc 
continuous  between  elements  while  the  pressure  is  approximated  by  a  linear 
function.  The  pressure  may  be  continuous  or  discontinous  between  elements.  The 
assembly  of  the  global  matrix  equations  is  by  the  direct  stiffness  approach  where  the 
equations  for  nodes  common  to  adjacent  elements  arc  added.  The  majority  of  the 
computational  time  is  spent  solving  the  nonlinear  equation  (43).  NACHOSII  uses  the 
Newton  method  to  linearize  the  equation  by  computing  the  Jacobian  of  the  system. 
More  details  on  the  implementation  of  this  method  will  be  presented  in  section  4.9.5. 

The  Newton  method  and  the  FDM  reduce  the  problem  to  a  linear  system  that 
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is  nonsymmctric  (for  the  nonlinear  Navicr-Stokes  problem)  and  indefinite.  The 
solution  method  in  NACHOSII  is  a  special  form  of  Gaussian  elimination  called  the 
frontal  solution  method.  This  method  combines  the  global  matrix  assembly  process 
with  Gaussian  elimination  so  that  unknowns  are  solved  in  the  system  as  soon  as  the 
assembly  for  each  element  is  completed.  This  combination  minimizes  the  amount  of 
memory  used  since  only  a  portion  of  the  global  system  is  ’active’  in  main  memory. 

For  large  problems,  e.  g.,  three-dimensional  problems,  iterative  linear  solvers  (e.  g. 
Conjugate-Gradient)  are  more  efficient  than  Gaussian  elimination.  Since  such 
iterative  systems  require  the  assembly  of  the  global  matrix  first,  no  attempt  was 
made  to  modify  the  frontal  assembly  process  to  allow  the  use  of  iterative  solvers. 

Besides  in-house  projects  at  Sandia,  NACHOSII  has  been  successfully  used  for 
the  solution  of  a  variety  of  problems.  NACHOSII  results  compared  favorably  with 
FIDAP  results  in  Gartling  (1990).  Koch  and  Yuen  (1985)  used  a  modified  version  of 
NACHOSI  to  model  mantle  convection. 

4.9.2  The  Inteerated  Method 

Recall  that  this  method  solves  the  continuity  equation  as  part  of  the  global 
system  of  equations.  Due  to  shortcomings  in  the  method  of  solving  the  nonlinear 
system  of  equations  (43),  this  method  did  not  keep  the  discrete  divergence  small 
enough  to  allow  an  accurate  solution.  The  problem  with  the  method  of  solution  will 
be  addressed  in  section  4.9.5.  Typically,  the  discrete  divergence  would  remain  small 
for  a  period  of  time,  then  increase  sharply  (usually  by  at  least  two  orders  of 
magnitude)  before  an  error  would  occur.  Figure  16  shows  the  change  in  discrete 
divergence  for  a  case  where  T^^,  =  25C  using  the  QUAD8  element  with 
discontinuous  pressure  approximation.  Following  the  large  increase  in  discrete 
divergence  at  step  269,  the  temperature  increased  to  an  unrealistically  large  value. 

This  variation  in  discrete  divergence  is  even  more  pronounced  in  the  T^^^  =  O.IC  case. 


82 


4.9.3  The  Penalty  Method 

Table  8.  How  e  affects  discrete  divergence. 


Penalty  Parameter 

Discrete  Divergence 

1.0  X  10*^® 

-10*® 

1.0  X  10**^ 

-10*® 

1.0  X  10*^® 

-10*® 

As  discussed  earlier,  the  penalty  method  uses  a  fixed  parameter  to 
approximate  the  incompressibility  condition.  The  choice  of  this  parameter  is  often 
difficult  since  choosing  it  too  small  will  ill-condition  the  problem  and  choosing  it  too 
large  will  introduce  too  much  compressibility  into  the  solution.  Table  8  illustrates 
how  changing  the  penalty  parameter  will  change  the  discrete  divergence.  Note  that 
these  calculations  were  for  the  T^^j,  =  25C  case  using  the  dimensional  form  of  the 
equations  and  a  QUAD8  element  with  discontinuous  pressure  approximation;  similar 
results  (at  smaller  discrete  divergences)  hold  for  a  properly  scaled  system.  The 
method  used  to  determine  the  appropriate  penalty  parameter  was  to  make  several 
sample  runs  and  observe  which  number  gave  the  best  discrete  divergence  behaviour. 
As  a  result,  early  steps  in  the  calculations  showed  great  promise,  but  the  discrete 
divergence  would  increase  sharply  at  some  point  later  in  the  calculation.  When  this 
increase  occured,  the  solution  soon  encountered  an  ill-conditioning  error.  As  a  result 
of  these  errors,  no  simulation  ran  to  the  point  of  showing  a  liftoff  of  the  ball. 

Figure  17  illustrates  this  performance  for  a  case  where  T^^j,  =  O.IC  and  €  *  1.0  X  10* 
using  a  QUAD8  element  with  discontinuous  presssure  approxmation.  The  penalty 


method,  however,  did  lower  computing  time  when  compared  to  the  mixed  method 
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since  it  solves  a  matrix  system  for  fewer  degrees  of  freedom.  The  penalty  method 
typically  used  about  35%  less  time  than  the  mixed  method  on  the  same  problem. 

4.9.4  The  Uzawa  Aleorithm 

The  Uzawa  algorithm  is  a  specific  application  of  an  augmented-Lagrangian 
method  to  solve  the  minimized  form  of  the  plume  problem.  An  augmented- 
Lagrangian  method  adds  an  additional  (or  augmented)  term  to  the  momentum 
equation.  The  advantage  of  this  method  is  that,  due  to  the  presence  of  this 
additional  term,  the  exact  solution  of  the  Stokes  problem  may  be  found  without 
making  the  penalty  term  go  to  zero  (with  the  resulting  ill-conditioning).  It  contains 
an  iterative  procedure  that  corrects  the  velocity  field  if  the  discrete  divergence  is 
too  large.  Detailed  theoretical  background  for  this  method,  with  applications  for 
Navier-Stokes  equations  may  be  found  in  Fortin  and  Glowinski  (1983).  Note  that 
they  mention  the  application  of  the  Uzawa  algorithm  to  the  time-dependent  Navier- 
Stokes  problem,  but  they  show  no  supporting  calculations.  In  fact,  no  application 
with  supporting  calculations  for  the  time-dependent  Navier-Stokes  problem  has  been 
found  in  the  literature. 

In  the  Uzawa  method,  the  matrix  equation  (38)  has  an  additional  term: 

M—*C{U)U-QP*KiU,f)U*B(,f)f*-Q^QU  =  F(7) 
dx  e 

where  the  1/c  Q^QU  term  is  the  augmented-Lagrangian  term  and  the  energy 
equation  is  unchanged,  e  is  a  small  constant  that  is  usually  called  a  penalty 
parameter.  This  method  has  the  same  advantage  as  the  penalty  method  in  that  it 
eliminates  pressure  as  an  unknown  in  the  matrix  system  and  thus  decouples  velocity 
and  pressure.  This  decoupling  reduces  the  size  of  the  matrix  system  and  saves 
memory  and  computing  time.  However,  the  pressure  solution  must  be  found 
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iteratively.  Now,  to  decouple  U  from  P,  start  with  an  initial  guess  for  P  and  solve 
for  U,  then  update  P  and  solve  again  if  the  discrete  divergence  is  too  large,  i.  e., 

1.  Given  1^"*,  solve  for  t)",  "f" 

M—  +  C(^")t7"  -  +  B(7*^7"  +  KO"'  =  6 

dt  ' 

N—  +  D(U*)f'  +  L(7*)f*  =  6 
dt 


where 


K,  = 

'  e  ' 

{  ff 

*  Jo 


2.  Update 


e  ' 

(This  step  is  based  on  a  perturbation  to  the  continuity  equation.) 

3.  If  discrete  divergence  is  too  large  repeat.  (The  traditional  penalty  method 
occurs  in  the  method  when  Pq  =  0  and  steps  2  and  3  are  performed  only  once.) 

Note  that  the  augmented-Lagrangian  theory  supports  this  choice  of  the 
additional  term  in  the  step  1  equation  only  for  the  linear  Stokes  problem,  i.  e.,  C  =  D 
=  0.  Since  the  linear  problem  can  be  expressed  as  an  equivalent  minimization 
problem,  the  minimization  theory  indicates  the  appropriate  choice  of  the  extra  term. 
The  nonlinear  momentum  and  energy  equations  that  describe  the  plume  problem 
have  no  known  equivalent  minimization  form.  Experience  has  shown,  however,  that 
using  the  same  additional  term  as  the  cv.-rresponding  Stokes  problem  produces  results 
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that  are  comparable  to  results  using  the  mixed  or  integrated  method  (Cuvelier,  et.  al. 
(1986),  Fortin  and  Glowinski  (1983),  Fortin  and  Fortin  (1985)).  This  Uzawa  method, 
incidentally,  is  not  difficult  to  implement  in  most  finite-element  programs,  since  all 
of  the  component  matrices  are  available. 

The  Uzawa  algorithm  seemed  very  promising,  particularly  in  light  of  the 
excellent  results  shown  by  Pelletier,  et.  al.  (1989).  They  stated  that  for  small  €, 
convergence  of  the  discrete  divergence  values  to  0(c)  occurred  with  only  one  or  two 
iterations  in  the  Stokes  flow  problem  that  they  solved.  Rapid  convergence  is 
necessary  for  this  method  since  the  global  system  must  be  solved  for  each  penalty 
iteration.  Since  this  solution  is  the  most  expensive  part  of  the  solution,  if  more  than 
two  or  three  iterations  are  used,  the  cost  advantages  of  the  penalty  method  are  lost. 
Table  9.  Sensitivity  of  Uzawa  to  changes  in  e 


€ 

Subiterations 
to  convergence 

Rate  of  Decrease  in 

Discrete  Divergence 

1.0  X  10-^^ 

I 

One  order  of  magnitude 
per  subiteration 

3.0  X  10*^® 

5-7 

18%  decrease  per 
subiteration 

1.0  X  10*^^ 

10 

10%  decrease  per 
subiteration 

Initial  testing  of  the  Uzawa  algorithm  showed  very  good  results  for  short 
periods  of  time.  The  value  of  c  greatly  affected  the  performance  of  the  algorithm. 
On  a  benchmark  plume  problem  with  a  15  X  20  grid  and  At  =  .0003  (50  steps)  for  the 
^amb  25C  parameters,  e  =  1.0  X  10'^^  required  ten  iterations  per  time  step  to  reduce 
the  discrete  divergence  to  0(c).  On  the  same  problem,  e  =  3.0  X  10*^®  required 
between  5  and  7  iterations  per  time  step  while  c  =  1.0  X  10'^®  required  only  one 
iteration.  This  sensitivity  is  also  reflected  in  the  rate  of  convergence.  When  the  c  = 
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1.0  X  10'^®  case  required  more  than  one  iteration  (after  100  time  steps),  only  two  or 
three  iterations  were  required  since  the  discrete  divergence  decreased  by  roughly  one 
order  of  magnitude  per  iteration.  In  contrast,  the  e  =  3.0  X  10"^*  case  saw  a  decrease 
in  discrete  divergence  of  about  18%  per  iteration  while  the  e  =  1.0  X  10*^^  case 
droppped  by  about  18%  per  iteration  .  These  results  are  summarized  in  table  9. 

These  observations  arc  consistant  with  results  reported  by  Pelletier,  et  al.  (1989)  and 
Fortin  and  Glowinski  (1983).  The  rapid  convergence  for  small  e  also  confirmed  that 
the  coding  modifications  for  the  Uzawa  algorithm  were  correct. 

Extensive  testing  of  the  Uzawa  algorithm  for  the  plume  problem,  however, 
consistently  demonstrated  that  the  performance  deteriorated  as  time  progressed. 
Simulations  that  started  with  ony  one  or  two  Uzawa  iterations  per  time  step  would 
often  required  4  or  5  after  100-200  time  steps.  A  typical  case  was  a  simulation  of  the 
^amb  "  25C  experiment  on  a  25  X  30  grid  with  At  =  .01.  The  first  70  time  steps 
usually  required  only  one  Uzawa  iteration;  only  on  3  occasions  were  3  iterations 
used.  The  subsequent  50  time  steps,  however,  required  5  or  6  iterations  per  time  step. 
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This  excessive  number  of  Uzawa  iterations  makes  this  method  less  efficient  than  the 
integrated  FEM.  For  this  case,  120  steps  using  the  Uzawa  algorithm  required  48 
minutes  of  cpu  time;  the  identical  case  ran  500  steps  using  the  integrated  method  in 
just  68  minutes.  Table  10  illustrates  how  the  rate  of  convergence  changes  as  the  time 
integration  progresses.  Note  that  these  results  are  for  another  sample  case  where  a 
15  by  20  grid  was  used  for  the  T^^^  =  O.lC  case  and  e  =  1.0  X  10'^*.  These  results  are 
presented  graphically  in  figure  18.  Since  the  convergence  was  highly  dependent  on 
the  choice  of  e,  and  a  single,  acceptable  value  for  e  was  never  determined.  These 
results  coincide  with  comments  made  by  Fortin  and  Glowinski  (1983)  for  similar 
nonlinear  problems.  Without  an  adequate  choice  of  e,  the  convergence  was  too  slow 
to  be  economical  when  compared  with  the  mixed  method.  These  poor  results  led  to 
abandonment  of  this  technique  as  described  here.  Other  similar  methods,  perhaps 
involving  an  automatic  choice  of  c,  may  be  more  effective  in  the  future. 

4.9.5  Multistep  Newton  Method 

NACHOSII  uses  Newton’s  method  to  solve  the  nonlinear  system  of  algebraic 
equations  that  are  formed  in  the  FEM  (equation  (43)).  Newton’s  method  gives 
quadratic  convergence  to  the  solution  provided  that  the  initial  guess  is  close  enough 
to  the  solution.  Since  NACHOSII  uses  a  predictor/corrector  type  time  integration 
(either  forward/backward  Euler  or  Adams-Bashforth/trapezoid  integration),  the 
predicted  solution  value  will  generally  be  close  to  the  correct  solution  for  each  time 
step.  This  assumption  provides  the  basis  for  NACHOSII  to  use  a  one-step  Newton 
method  for  the  solution  of  the  nonlinear  system.  This  method  is  justified,  with 
examples  given,  in  Gresho,  Lee,  and  Sani  (1980). 

The  modification  made  to  NACHOSII  is  a  simple  extension  of  the  base  logic. 
After  the  first  Newton  iteration,  the  discrete  divergence  is  used  as  a  test  of 
convergence.  If  the  discrete  divergence  is  too  large,  complete  another  Newton 
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Figure  18.  Discrete  divergence  for  Uzawa  method. 
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iteration  and  check  the  divergence  again.  Repeat  until  the  discrete  divergence  is 
small  enough. 

This  method  gave  acceptable  results,  as  reported  in  the  next  section.  Since 
most  time  steps  normally  required  only  one  Newton  iteration,  the  additional  cost  was 
not  significant  compared  to  the  original  one-step  Newton  method.  Typically,  those 
time  steps  requiring  additional  Newton  iterations  needed  only  two  or  three  iterations 
to  reduce  the  discrete  divergence  to  a  small  value.  A  plot  of  discrete  divergence 
versus  time  for  this  method  is  given  in  figure  19  for  a  15  X  50  grid  in  the  T^^  = 

O.IC  case.  Table  11  illustrates  the  rate  of  convergence  for  this  method  when  multiple 
subiterations  per  time  step  arc  required. 

4.10  Aoplication  to  Plume  Problem 

The  arguments  and  computational  data  presented  above  suggested  that  the 
Table  11.  Typical  convergence  for  multistep  Newton  method. 


Subiteration 

Discrete  Divergence 

1 

3.31  X  10'® 

2 

3.81  X  10-^^ 

best  approach  available  for  solving  the  plume  problem  was  to  use  a  QUADS  element 
with  discontinuous  pressure  approximation  in  a  mixed  FEM  that  employed  a 
multistep  Newton  method  to  solve  the  nonlinear  matrix  problem.  The  gjal  of  this 
numerical  work  was  to  develop  a  model  that  would  predict  the  laboratory  plume 
behavior.  The  following  results  will  show  that  this  method  ga\c  a  very  accurate 
prediction  of  the  laboratory  results. 

4.10.1  Ambient  Temperature  =  25C 

The  mesh  used  for  the  simulation  was  a  15  X  50  (figure  20a).  The  node 
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spacing  in  the  radial  direction  was  most  dense  near  the  axis  of  the  heater;  the  node 
spacing  in  the  height  direction  was  even.  For  all  cases,  the  maximum  element  size 
and  time  step  was  determined  by  the  fundamental  length  and  time  scales  used  or  the 
problem.  Based  on  the  theoretical  accuracy  of  the  FEM,  these  choices  also  insured 
that  the  numerical  resolution  was  greater  than  the  error  in  the  experimental 
measurements.  The  heater  was  represented  by  specifying  the  heat  flux  along  the 
bottom  of  the  first  four  elements.  Euler  time  integration  (forward  Euler 
predictor/backward  Euler  corrector)  was  used  early  in  the  problem  since  its  damping 
characteristics  removed  temporal  oscillations  in  the  solution  during  the  rapid  heater 
initiation  (the  trapezoid  method,  in  contrast,  has  no  damping  characteristics). 

Smaller  time  steps  were  used  during  this  time  due  to  the  reduced  accuracy  of  Euler 
integration  compared  to  trapezoidal  integration.  Later  in  the  problem,  larger  time 
steps  were  used  when  the  second-order  time  integration  method  (Adams-Bashforth 
predictor/trapezoid  corrector)  was  used. 

The  best  tool  for  plume  visualization  is  temperature  contours,  or  isothermal 
lines.  Some  typical  plots  are  given  in  figure  20b.  These  plots  also  provide  the  most 
effective  way  to  compare  the  numerical  results  with  the  experimental  results.  Figure 
20c  shows  the  results  of  plotting  the  height  versus  time  for  the  top  of  the  ball  for  the 
experimental  data  (top  plot)  and  the  numerical  results  (bottom  plot).  For  this  plot, 
the  top  of  the  plume  was  defined  by  the  T  =  25C  contour.  The  numerical  plot  is 
shifted  later  in  time  when  compared  to  the  experimental  plot  due  to  the  slower 
heater  start-up.  Note  that  both  plots  display  a  relatively  constant  speed  throughout 
most  of  the  flight.  For  the  experimental  data,  this  average  speed  is  .055  cm/sec;  for 
the  numerical  simulation,  this  average  speed  is  .056  cm/sec.  These  results  agree 
within  the  error  tolerance  of  the  experimental  data. 

4,10.2  Ambient  Temperature  =  O.lC 


93 


A  15  X  50  mesh  similar  to  the  T^jj«25C  case  was  used  for  this  ease  (figure 
21a).  Note  that  the  grid  is  placed  on  the  same  physical  dimensions  as  the  *  25C 
grid.  The  coordinate  values  are  different  due  to  the  change  in  the  length  scaling. 

The  element  spacing  in  the  radial  direction  is  different  from  the  T^j^  =  25C  case. 
This  spacing  was  chosen  to  place  smaller  elements  near  the  border  of  the  larger  ball. 
The  same  time-integration  approach  was  used  for  this  case  also.  The  only  difference 
between  this  case  and  the  previous  case  was  in  the  heater  start-up.  A  time-dependent 
heat  flux  boundary  condition  was  used  along  the  bottom  of  the  first  three  elements. 
The  heat  flux  was  increased  gradually  over  a  time  period  of  1.3  hours  to  the  full 
power  used  in  the  experiment.  This  gradual  increase,  along  with  the  damping  from 
the  Euler  time  integration,  smoothed  large  temporal  oscillations  in  the  early  steps  of 
the  solution.  This  slow  heating  also  shifted  the  plots  of  height  versus  time  further 
apart  than  in  the  case  (figure  21c).  For  this  plot  the  top  of  the  plume  was 

defined  by  the  T  =  15.4C  contour.  The  average  speed  results,  however,  agree  to 
within  experimental  accuracy-the  speed  of  the  experiment  was  .0055  cm/sec  and  the 
numerical  results  gave  .0053  cm/sec.  Figure  21b  shows  typical  temperature  contour 
plots. 

4.10.3  Ambient  Temperature  =  -26. 1C 

No  successful  runs  were  completed  for  this  case.  The  difference  between  the 
size  of  the  terms  of  the  global  matrix  corresponding  to  viscous  forces  (matrix  K  in 
equation  (42))  created  too  much  round-off  error  and  let  to  ill-conditioning  errors  in 
the  solution.  This  difference  is  not  due  to  scaling  as  discussed  earlier,  but  due  to  the 
wide  variance  in  the  value  of  viscosity  between  fluid  inside  and  outside  the  plume. 
The  difference  was  as  high  as  10^  for  this  case. 

A  method  that  may  have  promise  for  resolving  this  problem  in  the  future  has 
been  suggested  by  Pelletier  et.  al.  (1989).  They  suggested  scaling  the  matrix 
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equations  on  an  element  level.  This  scaling  could  show  a  similar  improvement  as  the 
global  scaling  of  the  full  equations  did  in  the  earlier  discussion. 
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Experimental  Result 
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Figure  20c.  Results  for  case. 


Figure  21b.  Temperature  contours  for  case,  time*3.1  hours. 
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Experimental  Result 
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Numerical  Result 
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Figure  21c.  Results  for  T^^j,=0.lC  case. 


CHAPTER  5 

SUMMARY  AND  CONCLUSIONS 

5.1  Experimental  Study 

The  experiment  described  in  chapter  2  explored  the  interaction  of  buoyancy 
and  diffusion  at  higher  viscosity  contrasts  than  earlier  studies.  Although  the  method 
of  plume  generation  (heater)  differed  from  the  related  experiment  (injection  of 
warmer  fluid)  of  Griffiths  and  Campbell  (1990),  the  plume  structure  appears  to  be 
dynamically  similar  to  their  results.  The  experiment,  then,  may  be  viewed  as  an 
extension  of  their  work  to  larger  viscosity  contrasts. 

The  results  of  the  experiment  (figures  4a,  5a,  and  table  3)  demonstrated  that 
previous  analytical  models  are  inadequate  for  describing  the  motion  of  the  ball.  The 
flow  structure  of  the  T^j,  =  -26.1C  case  (figure  6b)  suggested  that  heat  loss  from  the 
ball  is  an  important  factor  that  may  account  for  some  of  the  discrepancies  in  the 
earlier  models. 

Two  improvements  in  the  experimental  apparatus  and  procedure  would  make 
the  study  more  effective.  First,  a  more  powerful  heater  would  provide  the  ability  to 
study  the  effects  of  heater  power  on  the  plume  structure  and  speed.  The  heater  used 
in  the  experiment  of  chapter  2  required  full  power  to  produce  the  tracer  bubbles. 
Second,  accurate  temperature  measurements  would  improve  the  utility  of  the  results 
markedly.  Measurement  of  the  temperature  of  the  heater  surface  would  indicate 
what  rate  of  heater  startup  should  be  used  in  the  numerical  model.  The  temperature 
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on  the  ball  perimeter  would  indicate  which  temperature  contour  in  the  numerical 
model  to  use  to  mark  the  plume  boundary.  These  improvements  would  allow  for 
more  accurate  comparison  of  the  experimental  results  with  the  numerical  model. 

5.2  Analytical  Modelling 

The  efforts  in  analytical  modelling  focused  on  refining  the  Griffiths  and 
Campbell  (1990)  model  by  adding  the  effects  of  thermal  detrainment.  The  relative 
amount  of  mass  entrainment  and  heat  loss  was  varied  by  introducing  a  heat  loss  term 
that  depends  on  a  parameter  p.  For  all  finite  values  of  this  parameter,  the  new 
model  is  more  effective  in  modelling  the  experiment  than  previous  models. 

The  next  steps  in  improving  the  analytical  model  will  involve  relaxing  the 
assumption  that  heat  loss  in  the  conduit  is  negligible  and  refining  the  relationship 
between  the  entrained  fluid  layer  and  the  thermal  boundary  layer.  Calculations  of 
heat  flux  from  the  numerical  models  suggest  that  the  heat  loss  in  the  conduit  may  be 
as  large  as  30%  of  the  heat  loss  from  the  ball.  This  percentage  will  increase  as  the 
conduit  length  grows.  Adding  this  effect  to  the  energy  equation  in  chapter  2  will 
improve  the  accuracy  of  the  model,  particularly  for  longer  flight  times.  Recall  that 
to  develop  the  relationship  in  equation  (16),  we  assumed  a  linear  temperature  profile 
through  the  thermal  boundary  layer.  Making  this  profile  more  realistic  will  result  in 
algebraic  complication,  but  may  also  improve  the  model. 

5.3  Numerical  Modelling 

By  reducing  the  round-off  error  and  using  the  discrete-divergence  as  a 
diagnostic  to  choose  the  best  element/interpolation  combination,  numerical  results 
that  agreed  with  the  experimental  data  to  within  the  experimental  errors  were  found 
for  two  experimental  cases  by  making  a  small  modification  to  the  nonlinear  matrix 


103 


solution  procedure.  In  the  ■  O.IC  case,  the  viscosity  contrast  of  10®  was  higher 
than  in  previous  numerical  studies  of  plumes.  These  improvement,  however,  did  not 
allow  successful  simulation  of  the  -  -26.1C  case  due  to  the  large  differences  in 
the  magnitude  of  elements  of  the  diffusion  matrix  K  because  of  the  large  viscosity 
differences. 

The  successful  simulations  in  the  T^^  =  25C  case  and  the  T^i,  =  O.IC  case 
will  serve  a  benchmark  runs  for  future  improvements  in  the  method  of  numerical 
solution.  Improvement  efforts  will  concentrate  on  modelling  larger  viscosity 
contrast  and  in  making  the  solution  more  efficient.  The  use  of  element  matrix 
scaling  may  reduce  the  large  discrepancies  in  the  size  of  the  elements  of  the 
diffusion  matrix.  A  more  sophisticated  augmented-Lagrangian  method  may  reduce 
the  computational  time  required  for  the  problem.  Finally,  since  we  know  roughly 
where  the  large  temperature  gradients  arc  as  the  plume  develops,  some  method  of 
moving  a  fine  grid  with  the  ball  and  neck  would  reduce  computational  time 
significantly. 

An  improved  numerical  model  will  allow  more  detailed  applications  of 
numerical  experiments  in  wider  parameter  ranges  than  the  laboratory  allowed.  For 
example,  running  a  numerical  simulation  in  a  taller  tank  will  help  determine  if  the 
new  analytical  model  continues  to  work  well  for  long  flight  paths.  By  developing 
better  flow  visualization  tools  for  this  problem  (e.  g.,  some  way  of  indicating  streak 
lines),  we  should  be  able  to  find  a  way  to  measure  the  relationship  between  the 
thermal  and  fluid  boundary  layers,  and  thus  determine  P  from  numerical 
experiments.  The  overall  goal  of  the  future  application  of  an  improved  numerical 
model  is  to  simulate  mantle  plumes  and  thus  test  the  analytical  results  in  this  case. 
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FINITE-ELEMENT  EQUATIONS  FOR  THE  THERMAL  PLUME  PROBLEM 


Substitute  the  interpolation  expressions  for  velocity,  pressure  and 
temperature  into  the  governing  equations  (32)  -  (36). 

Momentum: 
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Continuity: 


*  0  (A2) 

' 


Energy: 
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ae  jl 


_a 

dx. 


The  Galerkin  FEM  minimizes  the  residuals  by  making  them  orthogonal  to  the 

interpolation  functions.  Take  the  inner  product  of  equations  (A1)-(A3)  with  the 
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interpolation  functions  to  get 
Momentum: 


r  l^“  a.,  r 
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Continuity: 


Q 


5 


(A5) 


Energy: 


(A6) 


In  arriving  at  equations  (A4)-(A6),  we  used  the  divergence  theorem  to  reduce  the 
second-order  diffusion  terms  in  (Al)  to  first-order  terms  plus  a  boundary  integral. 

We  also  allow  the  variable  material  properties  in  (A4)-(A6)  to  have  explicit 
spatial  variation  within  an  element  as  follows: 
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n(x^)  =  f  ^■(x^AW 
Kx^)  -  V(x)lait) 
a(x^)  X  V^(x)a(t) 

where  7  is  an  interpolation  function  and  the  ^  quantities  are  vectors  of  nodal  point 
material  properties.  Since  the  dependent  variables  are  known  at  the  nodes,  this  is  an 
effective  means  of  evaluating  material  property  values. 

The  integrals  in  (A4)-(A6)  may  be  evaluated  via  numerical  quadrature  to 
produce  coefficient  matrices  for  the  following  system: 

Momentum  and  Continuity: 
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Energy: 
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For  this  problem,  if  wc  separate  the  buoyancy  expression  from  the  general  force 
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vector  in  F,  we  can  write  (A7)  and  (A8)  as 


CiU)U  -  QP  *  K(G,f)U  *  B(f)f  =  F 


-qW  =  6 
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and 
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